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SUMMARY 


This study comprises two separate tasks: (1) an attempt to identify and 
minimize the uncertainties and potential inaccuracies of the NASA Multilayer 
Diffusion Model (MDM) using data from selected Titan III launches, and (2) a 
systematic analysis of the physical/chemical processes which take place during 
the buoyant rise of a rocket exhaust ground cloud and formation of a realistic 
time-dependent model. The former study is based on detailed parametric calcu- 
lations using the MDM code and a comparative study of several other diffusion 
models, the NASA measurements, and the MDM. 

The comparative studies and parametric calculations show that: (1) The MDM 
consistently overpredicts the ground level concentrations in the cases examined 
if the appropriate input standard deviation of the wind azimuth angle is chosen. 
(2) The current lack of micrometeorological information at KSC causes wide un- 
certainties in the results calculated from the model. (3) Environmental hazard 
analyses which require, as input, pollutant concentrations throughout the en- 
tire planetary boundary layer (PBL) must employ MDM with caution. (4) In a 
shallow PBL condition, the uncertainties in the entrainment constant used in 
the MDM cloud rise formula can cause a factor of three (* 3) error in the 
ground level predictions. For such conditions, the center of the exhaust cloud 
should be placed below the PBL height for the model calculations in order to 
give a conservatively high pollutant level prediction at the ground. (5) The 
other models included in the study contribute to an understanding of the rela- 
tive merits of the diffusion modeling technique. The strengths of these models 
can be used as guidelines for developing a new and advanced diffusion model in 
the future. 

The results of the second task can be summarized as: (1) The value 2500 

cal g -1 used in the MDM for the heat release from the Shuttle (or Titan III) 

solid fuel is reasonable. (2) Deluge water injected into the exhaust plume has 

little effect on the subsequent concentration predictions. (3) An average 

2 

loading of about 1 g m of alumina (about 15% of the alumina initially in the 
cloud) will be deposited on the ground in a worst case example, as for the 
Titan III launch of 20 May 1975. 



I , INTRODUCTION 


A, STATEMENT OF THE PROBLEM 

Concerns over the environmental impact of the toxic exhaust ground 
clouds generated by Space Shuttle launches have stimulated the development 
of mathematical models which can predict the dispersion of the exhaust 
effluents in the planetary boundary layer (PBL) under various meteorological 
conditions. One of these models (developed by NASA to provide real-time 
predicting capability) is the Multilayer Diffusion Model (MDM) described in 
Ref. 1. Some calculated results of ground cloud rise, cloud path, and HC1 
downwind concentration using this model for a number of Titan III 
launches have already been reported and analyzed (see e.g., Ref. 2). However, 
the assumptions incorporated into the MDM have led to concern that there might 
be large uncertainties associated with the predicted results. Critical 
analyses of the MDM model and its potential deficiencies can be found in 
Refs. 3-6. Since it is anticipated that the MDM will be used as the basis for 
assessing the environmental hazards of NASA launch vehicles, particularly the 
Space Shuttle, it was deemed necessary to: (1) evaluate the magnitude of the 

errors in calculations of concentrations and dosage fields using the NASA/ 

MSFC MDM under various meteorological conditions and (2) provide a method by 
which the present model can be improved so that realistic launch constraints 
can be developed based on the improved model. In order to answer these ques- 
tions, the chemical and physical processes occurring during the early period 
before diffusion becomes dominant must be considered because they may seri- 
ously alter the initial input required for a diffusion calculation and may 
cause environmental problems in the near field of the launch pad due to the 
dry and wet deposition of particles from the cloud. 

B, BACKGROUND 

The transport, evolution and atmospheric interaction processes of 
the exhaust cloud formed during rocket launches are most conveniently treated 
when separated into two stages: first, the buoyant force-dominated cloud 

rise phase and second, the atmospheric turbulence and ambient wind-dominated 
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diffusion process. When the cloud is initially formed at the launching 
pad, large quantities of debris from the ambient environment are entrained 
in it. The large size debris particles from the pad surroundings coagulate 
with the smaller aerosols of A1 2 0 3 , formed from the solid rocket propellant, 
and carry much of the A1 2 0 3 to the ground during the earlier periods of cloud 
rise. This phenomenon is known as dry deposition. As the cloud continues 
to rise and cool, condensation of HC1 and H 2 0 onto the remaining A1 2 0 3 and 
H 2 0 aerosol particles may occur. Therefore, when the acid aerosol particle 
grows due to coagulation and condensation and finally falls out, part of the 
HC1 will be removed from the exhaust cloud (wet deposition) . This dry and 
wet deposition, in addition to causing a possible environmental problem in 
the near field of the launch pad ( ^ 5 km) , also affects estimations of the 
HC1 inventory and the size and number of particles present at cloud stabiliza- 
tion. Errors in these composition estimates are particularly important since 
the data are used as input to diffusion models (e.g., NASA MDM) by which 
downwind ground level pollutant concentration is computed. In addition, a 
copious quantity of water is injected into the exhaust plume in the first 
few seconds in order to protect the launch pad structures from heat and to 
dampen and reduce acoustical energy feedback to the Shuttle. This deluge 
water may affect the chemical and physical processes during the cloud rise 
and may be especially significant in the prediction of the cloud stabiliza- 
tion height. 

Diffusion modeling of the rocket ground cloud in the PBL is further 
complicated by the lack of an adequate means to describe the microstructure 
of the PBL, as well as by the obvious existence of irregularities in initial 
cloud shape, uncertainties in pollutant concentration distributions, and non- 
uniformities in wind field and boundary conditions. In general, mathemat- 
ical modeling of the diffusion transport takes either of two approaches , 
i.e. , one obtains either an analytical or a numerical solution to the diffu- 
sion equation; in either case, for a practical problem, the turbulent terms 
are closed by the use of gradient transport theory. Analytical solutions to 
the diffusion equation can be obtained only for very limited conditions; for 
example, if the flow field is stationary and homogeneous, and the mean wind 
P r °file is uniform, the off-diagonal dif fusivities may be ignored and the 
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mean concentration distribution for the instantaneous point source in the free 
flow field can be expressed by the well-known Gaussian distribution obtained by 
solving the diffusion equation. Although numerical approaches provide the capa- 
bility of fitting more complex conditions, their applicability and accuracy are 
often limited by the large computer core storage on run time required. However, 
the ultimate accuracy of either approach still relies on the parameters chosen 
to describe the turbulence mechanism and the information on the advective wind 
in the flow field. Reviews of turbulent diffusion studies can be found in Refs. 
7 to 9. More detailed descriptions of several advanced diffusion modeling tech- 
niques for practical calculations are given in Section IV. A. 

The currently used NASA MDM program is based on the classic Gaussian 
model coupled with diffusion parameters from Cramer . 10 Two primary techniques 
are employed in the MDM, viz. , the unlayered first-order technique or the 
layered technique, each of which has two models. The two unlayered models are: 
(i) the cylindrical distribution, "model 1", employed when describing the com- 
plete exhaust plume in the mixing layer and (ii) the ellipsoidal distribution, 
"model 3", employed for the ground cloud when meteorological conditions allow 
use of the assumption of a quasi-homogeneous mixing layer. The two layered 
models are: (i) the static-plume, "model 2", and the multilayered distribution, 

"model 4". Model 2 is used for diffusive transport at altitudes of 3000 to 
8000 m where the turbulent transport mechanism can be ignored. In model 4, the 
surface mixing layer is divided into homogeneous (in a statistical sense) layers 
with a well distributed source. Models 3 and 4 assume that the top of the mix- 
ing layer totally reflects the effluents and that the ground surface has a 
specified absorption coefficient. As mentioned above, the Gaussian model 
results from a very ideal fluid field which rarely occurs in the PBL of the real 
world. Its applicability thus fully depends on empirical validation for each 
practical problem. In addition, the Cramer diffusion coefficients used in the 
MDM to model the physics of atmospheric turbulent transport are derived from a 
series of experiments which were carried out in the lower portion of the PBL. 
Therefore they need to be validated for the diffusion study of a rocket gener- 
ated cloud which is distributed throughout the entire PBL. 


4 



To provide a data base for refinement and verification of transport 
models, an extensive monitoring program was performed after a number of Titan 
III launches (e.g. , Ref. 11). During the diffusion process, data were 
taken with both airborne and ground-based sampling. Since the results 
calculated from MDM do not provide a time history of the pollutant concentra- 
tion at a given position, a direct and systematic comparison between the avail- 
able data, (especially the airborne data) by which refinements or verifications 
could be made, has never been conducted. 

C. BRIEF DESCRIPTION OF APPROACH 

This research program has comprised an effort to ascertain the 
aforementioned uncertainties and potential inaccuracies of the MDM, as well 
as a systematic attempt to model the physical /chemical processes occurring 
during the cloud rise period. The exhaust and early cloud properties such 
as the post-afterburning chemical composition, effective heat release, deluge 
water effects, and overall physical properties during cloud rise are discussed 
in Section II. This phase of che study is needed for predicting specific con- 
centration distributions in the stabilized cloud as well as for modeling the 
coagulation, sedimentation, and condensation processes which take place during 
cloud rise itself. 

Parametric calculations using the MDM code to identify and/or mini- 
mize error limits due to uncertainties in the input data have been partially 
discussed in previous studies. 3 ’ 6 In this report, the sensitive input quanti- 
ties identified in those studies have been used parametrically in calculations 
for different meteorological conditions. Our emphasis, however, is on the 
problems associated with the lack of information on turbulence in the PBL at 
KSC, which has not been discussed previously; Section III summarizes the results. 

With respect to the accuracy of predictions of cloud diffusion, this 
program has included a comparative study of five different diffusion modeling 
techniques including the MDM and comparisons with available measurement data. 

The objective of this phase of the study has been to address the following 
questions: (1) Do the measurements conducted by NASA provide a useful data 

base for the verification of diffusion models? (2) What would be the predicted 
results if other models were used in making an assessment, i.e., can they yield 
better results than MDM? How? and Why? (3) Can guidelines be developed for 
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improving the performance of MDM without substantially altering its present con- 
dition? (4) What kind of approach for modeling the atmospheric diffusion of a 
rocket generated cloud in the PBL is indicated, if a better assessment is 
required in the future? 

The modeling techniques examined in this study, in addition to the MDM 

are : 

(1) The Air Force/Vandenberg METS 12 model which is a Gaussian model with 
Turner-Pasquill stability classes for the determination of power law 
type dispersion parameters. 

(2) The NUS TREATS 13 model which uses a Gaussian assumption for horizontal 
distribution and an integrated moment scheme for determining disper- 
sion parameters. 

(3) The Lawrence Livermore ADPIC 14 code which is based on a particle-in- 
cell/pseudovelocity transport scheme in which only diagonal eddy 
dif fusivities are used. 

(4) The AeroChem DISF 15 code which is a diffusion model for shear flow 
obtained by solving the diffusion equation via a second-order eddy 
diffusivity closure scheme and the Lagrangian approach for statisti- 
cal quantities. 

This study concentrates on the clouds produced by Titan III rocket launches at 
Kennedy Space Center on 10 December 1974, 20 August 1975, and 14 March 1976. 

The comparative results and the interpretation of meteorological information 
are given in Section IV. The reasons for selecting these three launches are: 

(1) more complete experimental data are available on these launches than for 
others; (2) the Titan III rocket uses the same propellant as the Shuttle SRM 
(and is about 1/2 the mass flow), and (3) these launches cover a broad range of 
interesting meteorological conditions. To elaborate on the latter point, it 
should be noted that the December launch was a night launch in the winter time 
with the presence of a shallow stable stratified PBL; the August launch took 
place in a daytime sea-breeze condition at KSC in which a high inversion layer 
was present; and the March launch was accompanied by strong humidity, variable 
wind direction and a moderately weak inversion layer (see Section IV. B for more 
detailed discussion) . 
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A model for the coagulation of AI2O3 and debris particles has been 
developed to provide particle size distribution, particle densities, the mass 
fraction of alumina and the mass of sedimentation in each particle size class 
prior to the onset of condensation. Section V is devoted to the formulation of 
this model, its validation and the results of calculations. 

The authors are indebted to Scott Wagner, Joseph Mathis, Richard 
Bendura, Richard Gomberg, Gerald Gregory, and Gerald Pellett of NASA/Langley 
(Environmental Field Measurement Branch, Marine and Applications Technology 
Division) for their constructive discussions and advice during this study. 

The authors also wish to thank the people who ran their diffusion models for 
this comparative study, in particular Christine Sherman and Rolf Lange at 
Lawrence Livermore Laboratory, Martin Chandler at NUS Corporation, and Darryl 
Dargitz at Vandenberg Air Force Base. The helpful comments and assistance on 
the report from Dr. William Miller at AeroChem are gratefully acknowledged. 

We also thank James Mills who coded the CRAM program and ran the computer pro- 
grams for us at AeroChem. 
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II. EXHAUST AND CLOUD RISE PROPERTIES 


In our overall approach to modeling the physical chemical processes 
in the cloud rise regime, the exhaust and early cloud properties have been 
calculated through the following rather complex but straightforward procedure. 

A. EXHAUST COMPOSITION 

The post-afterburning exhaust composition is calculated based on 
the assumption of complete combustion. An extensive study 16 using exact 
values for nozzle exit plane compositions and conditions and a low altitude 
afterburning/mixing plume program (LAPP 17 code) has shown that this assumption 
is very closely approached. This is especially true for the major species 
which are the only ones that need be considered since the chemistry of the 
minor species is not important in the determination of overall cloud 
characteristics . 

The post-afterburning plume composition is shown in Table I. The 
values given include added air but only that required for combustion. This 
amounts to about 149% (by weight) of the rocket propellant for the solid 
engine (SRM) and 137% of the propellant weight for the liquid engine (orbitor) . 
In addition. Table II gives comparative values for the masses of three major 
post-afterburning exhaust products that result from calculations invoking the 
above assumption and those produced by extensive calculations using the LAPP 
code for Titan and Shuttle SRMs . The values given are all based on the pro- 
ducts of 100 g propellant. It should be noted that the value cited by 
Gomberg and Stewart 16 for this comparison has been derived by using the 
amount of alumina (A1 2 0 3 ) in the plume as a normalization factor with which 
to obtain the amounts of other gases not specified in their calculations. 

They show a very good agreement if we account for the error and uncertainty in- 
duced from the thermal-kinetic data, in the numerical scheme, of the exit plane 
conditions, etc. 

B. HEAT CONTENT AND EFFECTS OF DELUGE WATER 

With the complete burning assumption, the heat release from the 
propellant is found to be 2847 cal g -1 for the solid engine SRM and 3117 cal g -1 
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for the liquid orbitor engine. These values do not include the effect of the 
deluge water. The effective fuel heat content for the solid engine is essen- 
tially the same as the value calculated by Stewart and Grose 3 and is slightly 
higher than the value of 2500 cal g- 1 currently used in the NASA/MSFC MDM. 1 The 
value for the liquid engine is about a factor of 6 higher than the value used in 
the NASA/MSFC MDM (version 5) . 1 

The principal potential effects of the deluge water on the 
chemical/physical process during cloud rise are: (i) changes in chemical 

composition due to possible quenching of some of the high temperature after- 
burning reactions and (ii) reduction of the effective heat content in the 
cloud due to the vaporization of the deluge water. A study 18 of the former 
problem was recently conducted at NASA/Langley and concludes that changes in 
the concentrations of major species are negligible. In other words, the 
complete burning assumption used in the present study is still applicable 
even during the first few seconds of firing in which the deluge water is 
injected. In considering the effect of the deluge water on the effective 
heat content and subsequently on the prediction of the cloud stabilization 
height using the present cloud rise model, we have included all of the water 
to be used in the Acoustic Water Suppression System (AWSS) for the launching 
of the Space Shuttle in order to determine its maximum effect. In the present 
plans for the AWSS, water will be poured under the firing pad and into the 
plume trench at a rate of 6500 gallons per sec over about 8 sec after rocket 
ignition. This represents a mass flow rate of about 2.46 x 10 7 g sec -1 . If 
the stabilized cloud is assumed to contain 17 sec of exhaust, which is a 
reasonable value based on Titan ground cloud data (a detailed discussion of 
this point is given below), then the total amount of effective available heat, 
will be 

Q h = M s x h s x t f + M L x h-L x tf - h w x M„ x 8 = 4.27 x 10 11 (cal) (1) 

where M s , M^ and M„ are the mass flow rates (g sec -1 ) of the solid engine, 

liquid engine, and poured water, respectively, h g and h^ are the effective 

heat release (cal) per gram of propellant, h is the latent heat of water 

w 

vaporization at room temperature, and t^ is the firing time during which the 
cloud is formed (now = 17 sec) . 
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This value for the total heat content equals about 80% of the total 
heat content excluding the effect of deluge water. Since the predicted stabil- 
ization height, Z^, is proportional to the 1/4 power of the total heat content 
(i.e., Z m <* Q h ^ 4 ) in the presently used Briggs 19 formula in MDM, the effect of 
deluge water on the prediction of the stabilization height is less than 5% in 
the worst case. It can certainly be ignored. For the Titan III rocket, the 
mass flow rate of trench water (which is the only deluge water in use) is about 
1.4 x 10 6 g sec -1 but it is employed over only a 5 sec period at most after fir- 
ing. Therefore, the average heat available for the first 10 sec of launch of 
Titan III is about 2750 cal g -1 if the trench water is considered. This value 
is slightly higher than that used in the NASA/MFSC MDM; the differences, however, 
are apparently negligible for the reason mentioned above. 

C, PHYSICAL PROPERTIES IN THE CLOUD RISE REGIME 

Since the net effect on buoyant cloud rise of atmospheric entrain- 
ment due to turbulence is not well known, an exact calculation of the cloud 
shape, its concentration distribution, and other physical properties during 
cloud rise would require an expensive and time-consuming direct numerical 
s im ulation process, Such a calculation is presently impractical. There- 
fore, we have continued to use the Briggs formula, as in the NASA/MSFC MDM, 
for the estimation of cloud rise speed when observations are not available. 

The cloud shape is assumed spherical; the rest of the physical properties, 
especially those required in calculations of coagulation and condensation 
(or heterogeneous chemistry) during the cloud rise, are obtained from a model 
based on an energy balance, combined with the above assumptions for cloud 
volume and cloud rise speed (Briggs formula) . The primary physical proper- 
ties required, in addition to the volume and cloud rise speed, are the temper- 
perature and mass of air entrained. 

The following diagram illustrates the formulation of the energy 
balance during the cloud rises 
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where Q is the total effective heat content in the cloud, P, P , are the 
n amb 

instantaneous pressure in the cloud and ambient environment where the cloud 

locates, m . is the mass of air entrained due to turbulent mixing and V , is 
axr & obs 

the observed cloud volume at the following stage. 

When the buoyant cloud rises, it will cool mainly due to expansion 
and entrained cool air. An instantaneous pressure balance between the cloud 
and the ambient environment can be reasonably assumed. With this assumption, 
the temperature and the mass of air entrained can be calculated, using the 
perfect gas law, if the real cloud volume is known. The mathematical form 
for this formal process can be written 

t T 

\-J pdv ' s D ( m i f c p. dT ') 

J ie r J 

to 298 


+ 



to 


air 


m . 
air 


( T (t ' ) - T amb (t'))dt' + 1/2 M W 2 (t) (2) 
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and 


t 

9 

M(t) = J m a±r (t')dt’ + m ± (3) 

to 

where P is the set of post-afterburning chemical species, in which the ele- 
ments are identified by means of a sub-index "i", W is the cloud rise speed, 
and m, T, t, and Cp represent the mass, temperature, time, and specific heat 
constant, respectively. The subscript "air" refers to entrained ambient air, 
"amb" to ambient, and "obs" to observed. 

The terms on the left hand side of Eq. (2) give the total heat 
content, and the work done by adiabatic reversible expansion due to the cloud 
rising from higher to lower ambient pressure. The terms on the right hand 
side consist of the internal energy of each species in the cloud, the energy 
required to heat up the entrained cool ambient air and the kinetic energy of 
the cloud. The potential energy is negligible. The total mass M(t) at time 
t in the second equation can be obtained from the observed volume, tempera- 
ture and pressure by using the perfect gas law; i.e., 


M/V , ^ -v *T (t) »K 
obs(t) gas 


' P a„b (t) 


(4) 


where K is a gas constant, i.e., K = n x R, 
gas & gas 

where R = 8.2 x 10 -5 (atm»m 3 /mole»K ) and n is the number of moles in a unit 


mass which equals 



where MW^ is the molecular weight and is the mass fraction of species i 
in the cloud. 

The total heat content, Q^, used in Eq. (2) is calculated from 
Eq. (1) for the Shuttle. If the observation data for the cloud rise history 
are not available, the cloud rise speed W is obtained from Briggs' formula 
which gives 

W = 1/4 A 1 ^ S l/i [sin(t s l / 2 )][l - cos(ts 1 / 2 )r /4 (5) 

where s = ■£- , the Brunt Vaisala frequency 

T r <3Z 
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where p a j_ r and T„ are the density and temperature, respectively, of ambient 
air near the ground, y is the entrainment constant and is the vertical gra- 
dient of ambient virtual potential temperature. Also if an observed volume is 
not available, the spherical-cloud shape assumption gives the volume as follows 


V obs^ = (Y*Z(t) ) 3 •tt« - 


and 


Z(t) 




W(t)dt + Z c 


( 6 ) 

(7) 


where Z(t) is the cloud altitude at time t and Z 0 is the initial cloud center 
location. t^ in Eq. (7) represents the time which will give the initial 
cloud rise speed from Eq. (5). 

The second term in the left hand side of Eq. (2) can be written in 
algebraic form, i.e.. 


h 


- _ _1 v p 1/Y' pl-l/Y' , I vp 

yt w or 0 r + yt v 0 Po 


( 8 ) 


where y' is the specific heat ratio, i.e. 


y' = —2. 

T Cv 


Substituting. Eqs. (3) to (8) into Eq. (2) closes the energy balance equation. 
Then the historical temperature variation during the cloud rise as well as 
the mass of entrained air can be solved implicitly. A numerical program 
following the above procedure is given in Appendices A and B as a subroutine, 
CLDRISE, of the cloud rise aerosol model (CRAM) program. 

A set of results using this model for volume, temperature rise, 
speed and entrained air mass for three Titan III launches at KSC is shown in 
Fig. 1. There is insufficient observation information at present to make 
meaningful comparisons. However, these can be made whenever data become 
available . 
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III. PARAMETRIC STUDIES USING THE MDM CODE 


A parametric study of version 5 of the MDM program was carried out to 

(1) place bounds on the uncertainties in MDM predictions of exhaust effluent 
concentrations and (2) establish guidelines for the use of the MDM and/or pro- 
vide directions for its modifications. Such a study also provides good back- 
ground for comparative studies using the various diffusion models employed in 
this program. 

The present work concentrates on the MDM model 4 because it is the 
model presumed capable of accounting for the inequalities of the source cloud 
(at stabilization) and the inhomogeneous nature of the real world. Since each 
species in the cloud is treated as an inert gas in MDM, it is sufficient to 
consider only one species, namely HC1. The meteorological conditions examined 
were those occurring during the two Titan III launches, viz., 10 December 1974 
and 20 August 1975. Figures 2 and 3 show the temperature, wind speed, wind 

direction, and humidity for these two launches. The meteorological data are 
adapted from the most recent rawinsonde measurement taken before the launch and 
are used as the reference basis for the calculations. 

The original parametric studies with MDM were made by Stewart and 
Grose. 3 Their work mainly covered two specific meteorological conditions, 
viz., low level sea-breeze conditions and the fall fair weather regime. 

Table III, adapted from their report, summarizes the kinds of parameters and 
the range of their variations in this work. Their reasons for selecting these 
parameters and this range of variations are not given here. Briefly, they con- 
cluded that the ground level concentration predictions using MDM are strongly 
dependent on the cloud geometry assumed at stabilization, the concentration 
distribution, the depth of the PBL, and the layer transition time. It should 
be noted that (1) their study was based on version 2 of the MDM 20 program, 

(2) the meteorological conditions they considered are different from ours, 
and (3) the operational procedure for preparing the input for MDM has been 
changed since then. It should also be noted that the initial work by Stewart 
and Grose did not cover all the parametric uncertainties in the MDM. Hence a 
further discussion, to a sufficient extent to support our later comparative 
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study, should be given. This report emphasizes parameters which have not 
been discussed in those previous studies. They are separated into two cate- 
gories; the first concerns those parameters which are related to the diffusion 
process and the second comprises those which affect the initial description 
of the cloud. A brief discussion of the calculational results using the sen- 
sitive input parameters found in previous studies is given in Section 2II.C., 


A. PARAME TERS AFFECTING DIFFUSION 

The thred kinds of parameters in this first category are the diffu- 
sion parameters, the boundary absorption factor, and the mean wind speed. 


1. Diffusion Parameters 


As mentioned above, the MDM program assumes a steady-state Gaussian 

dispersion. Using the Gaussian distribution form, the only unknown parameters 

are the standard deviations of the mean concentration dispersions, a , a , and 

x y 

o z , in alongwind, crosswind and vertical directions, respectively. The along- 
wind dispersion parameter o x used in MDM is based on the expression for along- 
wind cloud growth derived by Tyldesley and Wallington, 2 1 i.e.. 


a 


x 




2 


where is alongwind standard deviation of initial pollutant distribution 


and 


iL = 0.28 • ~~ • x 


A* 

as U > 0 or — - < 0 

Az 


L = 0 


for all other cases. 


where AU is the wind speed difference between that at the top and the bottom of 
the PBL. The problems associated with the use of this expression have been 
discussed by Stewart and Grose. 3 Since Dumbauld (one of the authors of the MDM) 
recommends using 0.6 rather than 0.28 in the above expression because in his 
experience, it yields more realistic results, we have used 0.6 in all calcula- 
tions in the present study. The resulting peak in the concentration maxima ob- 
tained with 0.28 is about 80% higher than that obtained with 0.6. 
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The expressions for crosswind and vertical dispersion parameters, 0y 

and a y in MDM are based on Cramer's formula, 10 in which they are related to 

the standard deviation of the wind azimuth angle a^(z,x) and elevation angle 

0g( z,t), where z is the height of source location and T is the so-called source 

release time. At present, t is taken as the time required for the cloud to rise 

to the stabilization height and z is the height of the center of each subcloud. 

In the preprocessor program, a (z,x) and a (z,x) in the entire PBL are 

A E 

determined by a very simple expression, i.e., 

a A (z,x) = a E (z,x) = a k J 2 (9) 

where a is the standard deviation of azimuth angle measured over 10 minutes at 
Ao 

a reference height (usually at 4 m above the ground) . 

Expression (9) implicitly incorporates the assumptions (1) that the 

atmospheric flow is isotropic, i.e., a = a and (2) that a and a are indepen- 

dent of height. These assumptions are inadequate because the atmospheric flow 

field is nonisotropic and o and are only independent of height near ground 

level (the altitude is unlikely to be more than 15% of the PBL thickness, see 

Section IV.B.l); in general, a and <J decrease monotonically at upper portions 

of the convection layer. It should also be noted that the assumption, = 

O / 2 in (9) has not been proven theoretically or empirically. In addition, 

Ao 

while direct measurements of c. have been made at each tower at KSC, the values 

A 0 

of a. vary considerably for different towers at different times. Therefore, 

Ao 

the effects due to the inadequate assumptions and the uncertainties of the value 

of O. should be determined. 

A 0 

Two approaches are utilized to determine the effects caused by the 
above problems. 

(i) Assume 

°a ct f 

CT A (z ’ t) = 2 and a E (z ’ T) = 2 

where a is the standard deviation of wind elevation 
Eo 

angle measured over 10 minutes at a reference height a 

few meters from the ground. The values of a and a 

Ao Eo 

are obtained from four sources: (1) the average value of 
the tower measurements, (2) the interpolation method 
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suggested by Record et al, 22 (3) similarity theory, see 

Section W.B.l, and (4) the value given by Susko and 

Stephens. 33 Table IV lists the values G f a an d a 

A 0 E 0 

obtained from these four sources. Since only a. 

A.Q 

is available for tower measurement , a simple assumption 

that a. = CL. is used in the calculation for source (1). 

Ao iIjo 

(ii) By definition, cr A (z,x) - O v (z,x)/U(z,x) and 

a (z,x) - cr (z,t) /U(z,t) , where a and a are the 

w V w 

standard deviations of the turbulent wind fluctuation 
components in the crosswind and vertical directions, 
respectively. Values of a ^ and are obtained from the 
method discussed in Section TV. B.l. Since the values of 
and ct e used in the MDM are dependent on the release 
time, we have used the expression given by Cramer et al, Zi * 
i. e. , 

a A (z,x) = a A (z,Xo) 

The predicted maximum ground concentrations from the five different 

sets of a, and o_ described in (i) and (ii) above for the December launch are 
A E 

shown in Fig. 4. As shown, they differ significantly in the near field 
(2 km i x S 10 km) of the launch pad. In particular, the values obtained 
using either the measurement data (i.l curve) or the theoretically more exact 
interpolation method for the a A and a E (ii curve) are much lower than the 
values obtained using other adopted methods. It should be mentioned that the 
predicted results from MDM could be even lower than the lowest value shown in 
Fig. 4 if one used the measurement data for a A and cr^ from the lower portion 
of the PBL and the monotonic decay profile for the upper portion. This is 
because (1) a vertically uniform profile for a A and is implicit in the 
calculation using measurement data, and (2) the a A and for the lower portion 
(calculation (ii) above) were about 7° and 4°, respectively, which are higher 
than the 2° measured value (see Tahle IV). However, it is clear that the 
present lack of turbulence information can lead to a wide uncertainty in the 
MDM predicted results. 
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2 . 


Earth Absorption Factor 


The MDM program (models 3 and 4) utilizes the simple mirror reflection 
method to treat boundary effects. An absorption factor, Yp is used to determine 
the strength of the absorption of the earth surface. For instance, water sur- 
faces may absorb all HC1 reaching them, but some land surfaces, like highways, 
may reflect most or all of the HC1. Therefore, we cannot classify the earth sur- 
face in the KSC area as being uniform. Since Yp is assumed uniform, it is of 
interest to see the effect of varying the value of Yp on the ground concentra- 
tion prediction. The results obtained from the two launches are almost identi- 
cal, as shown in Fig. 5. It can be expected that a greater effect would be 
observed in the far field, but even there the variation is less than a factor of 
two. Since the greatest variation in the peak of maximum concentration predic- 
tions (or 10 minute time-mean-concentration) is less than 50%, we can conclude 
that the surface absorption parameter is not an important factor in predicting 
ground pollutant level. For reference, comparisons of the results of varying 

Yr, and a. and a in MDM for the December launch are shown in Fig. 6; it is 

y Ao i^o 

obvious that the diffusion parameter dominates. 

3 . Wind Speed 

The wind profiles used in the diffusion calculation at KSC can be ob- 
tained from one of three measurements, viz., rawindsonde, windsonde or jimsphere. 
The rawindsonde measurements are used more widely probably because they provide 
more complete mesoscale meteorological information, (i.e., temperature, pressure, 
wind speed, wind direction, humidity, etc.,). However, all the measurements 
are taken at different locations and at different times, often not at the time 
of launch. Furthermore, even the wind profiles adopted from these three differ- 
ent measurements at very close intervals show considerable differences (see Fig. 
7). In other words, the wind profile used in the diffusion predictions may not 
represent the real situation and the effect caused by the uncertainty of wind 
speed must be determined. In the present study, calculations were made by 
changing the magnitude of wind speed but not the wind profile. It was found 
that in both the December and August cases no more than ± 85% change in the 
time-mean-concentration was produced by + 100% change in average wind speed, but 
there were no changes in the maximum concentrations. This is because the 
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formulation for the maximum concentration in the MDM program is independent of 

the magnitude of wind speed and the turbulence parameters , O . . and a , in the 

Ao Eo 


calculations were invariant. 


Theoretically, the a and a are dependent on 

a Ao tjo 

w 


wind speed: because o. — a. — and a and a are not dependent on the 

AU’AU’ v w 


magnitude of U, and Cg will vary with U. Since, in our calculations the a ^ 

or a are based on the tower data, there is no reason to change them when we 

£.0 

vary the wind speed obtained from soundings. 


B, PARAMETERS AFFECTING INITIAL CLOUD DESCRIPTION 

The second category of parameter studied includes the 

vertical mean gradient of ambient virtual potential temperature, A<J>/Az, the 

entrainment constant, y, and the standard deviations a , a and a for the 

Xo y 6 z o 

source distribution in the stabilized cloud. The effects of varying a are 

z o 

covered in the discussion on rearranging the vertical source strength. 

Section III. C. ; the remaining parameters will now be addressed. 


1. The Gradient of Virtual Potential Temperature Atfr/Az 

If A<j>/Az has a factor of two uncertainty, then a 20% variation 

for the stabilization height, Z^, can be induced, because “ 

For the present MDM program, a 20% change in stabilization height can cause 
"at most" a 230% increase (or decrease) in the total amount of pollutant con- 
sidered in the dispersion. Since the MDM program is an inert diffusion model, 
an increase (or decrease) of 230% in the initial source strength would cause 
a 230% increase (or decrease) in the predicted ground level concentration as 

well. Therefore a factor of two variation in A(j>/Az implies "at most" a 230% 
variation in the prediction of the ground level concentrations. However, 
since a big change in A<f>/Az would mean significant change in the atmospheric 
structure and thus question the validity of MDM, it is logical to assume that 
A<J>/Az is approximately constant over time and space. Also, the three factors 
for determining A<f>/Az, namely temperature, pressure, and humidity, are fairly 
reliable measurements and it is unlikely that Acf>/Az would have a factor of 
two variation. 


m\ ■ l/ “ 

Az J 
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2 . 


The Entrainment Constant y 


For the rocket generated cloud, the entrainment constant used in the 
Briggs formula in the MDM model has not been fully validated. However a pre- 
vious study by Hart, 23 using a one-dimensional cloud rise model* showed that 
the proper entrainment constant is about one on the basis of comparisons with 
the limited observed data of Titan III clouds. This entrainment constant value 
is different from the recommended value of 0.64 in the MDM model. It should 
also be noted that Hart's study showed very limited success in comparisons with 
the Titan III clouds in stable conditions. Therefore the effects due to the 
uncertainty of the entrainment constant should be investigated. Since, for 
stack plumes, the value taken is usually between 0.5 and 1.3, 2 6 we have made a 
series of calculations varying y from 0.5 to 1.3 for the December and August 
launches. The resulting ratio of source strength is given in Fig. 8. As 
shown, y has an insignificant effect on the diffusion prediction in the August 
case which is a daytime launch and we can generalize that the sensitivity of 

MDM to the entrainment constant in a deep PBL condition is fairly low. By con- 

x 

trast, the entrainment constant can cause a factor of three ( ± 3) variation in 
predicted ground level results for a shallow PBL. 

3. The Standard Deviations, Ao. and gy n , of Initial Pollutant 
Distribution 

The a Xo and 0y o now in MDM are chosen arbitrarily. In order 
to determine their sensitivity in the diffusion calculation, they have been 
varied consistently for all sublayers, by ± 10%, ± 30%, and ± 50%; such con- 
sistent variation keeps the elliptical cloud shape assumption, which is incor- 
porated in MDM, intact. The resulting peak of maximum concentrations for the 
December launch is shown in Table V as an example. It is found that there 
are insignificant changes in the predictions resulting from the above vari- 
ations in these two parameters for the December launch. However, there is 
about an 85% increase for the August launch when one makes a 30% decrease in 


The main difference between this model and the Briggs cloud rise formula is 
that this model introduces an additional equation of state to more accurately 
account for the buoyancy force; this equation of state is basically similar to 
the one described in Section II. B. 
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the a Xo and CTy 0 ; in contrast there is only a 35% decrease for a 30% increase of 
°x 0 and CTy 0 . It should be kept in mind that the sensitivity of the standard 
deviations of the initial pollutant distribution to MDM is higher for the case 
of a deep PBL than for a shallow PBL. 

C, SENSITIVE PARAMETERS IN PREVIOUS STUDIES 

Calculations using the sensitive parameters found in previous studies 
for the August and December launch cases have been performed. The resulting 
conclusion is similar to that of previous studies except that varying the source 

strength distribution in each layer would not have a significant effect on the 

ground level concentration if a notable amount of the exhaust were not mandator- 
ily placed in the lower portion of the PBL. In addition, it should be men- 
tioned that over a realistic range of variations in the depth of the PBL the 
change in the results is within a factor of two. As a precise example, although 

the depth of the PBL for the August launch cannot be determined exactly, it 

should range between 1000 and 2000 m. (Further discussion is given in Section 
IV. B. 2.) Under such consideration, the calculations showed that the variation 
of the predicted peak of maximum ground level concentration is less that a fac- 
tor of two. Nevertheless, since other studies 27 have shown that the peak in 
mass concentration at the ground level can reach 4 ppm while 8 ppm is the short 
term acceptable limit, the error of a factor of two caused by the uncertainty of 
the depth of the PBL cannot be neglected. A model which is capable of determin- 
ing the depth of the PBL is apparently desired. 
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IV. COMPARATIVE STUDY OF GROUND CLOUD DIFFUSION MODELS 


A. A BRIEF REVIEW OF SELECTED DIFFUSION MODELS 

Basically every model of atmospheric diffusion processes starts 
with the diffusion equation (formulated from the mass conservation equation) 
and then uses some technique to close the turbulent fluctuation terms in the 
equation. The usual technique, called K theory or gradient transport assump- 
tions, employs eddy dif fusivities to relate the gradient of mean quantities 

to the turbulent flux, e.g. , u c r = K,-j - 5 —. Here, u.’ is the fluctuating 

x dxj 1 

velocity component in the i direction, c' is the fluctuating concentration 
and the overbar denotes an average quantity. K-^j is the eddy diffusivity 
tensor which is obtained by utilizing either one or another of many possible 
assumptions or empirical results. Since the K theory has been virtually the 
only kind used to date for practical calculations of atmospheric diffusion 
processes, all modeling techniques selected for this comparison are limited 
to this category. Table VI lists these selected models, each of which is 
briefly described below. 


1. NASA/MSFC Multilayer Diffusion Models (MDM ) 

The MDM 1 program is based on the classic Gaussian distribution assump- 
tion (see Section I.B). The standard deviations of the mean concentration dis- 
persions, O x , Oy, o z , (usually called dispersion parameters) are obtained from 
Cramer's formula for vertical and crosswind terms and from Tyldesley and 
Wallington's 2 1 f or the alongwind term. The treatment of vertical wind shear, 
turbulent nonhomogeneity and source irregularity involves dividing the flow 
field (the PBL) into sublayers and having the source in each sublayer dis- 
sperse in a Gaussian form which depends on its initial local properties. A 
more detailed discussion of this model can be found in Ref. 5. 


2 . Meteorological Effluent Transport Simulation Model (METS) 

Using an approach similar to the MDM program, the METS model 12 em- 
ploys the Gaussian distribution assumption for the gaseous and liquid constit- 
uents and the layered structure. The vertical and the crosswind dispersion 
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parameters are given by power law expressions and the alongwind parameters 


equations , 

i.e.. 
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where a , CT are the initial standard deviations, and x is the downwind 
z 0 yo 

distance from the center of the source at x 0 . The exponents py, p z are deter- 
mined by fitting the empirical Turner-Pasquill 2 8 curves for different stability 
classes. The stability class for each sublayer is dependent on the calculated 
local gradient Richardson number (see Section IV. B). 

Two additional models have been incorporated into the METS model to 
account for the diffusion of particulates (e.g., A1 2 0 3 ) and the hydrochemi- 
sorption of HC1. These additions are beyond the scope of the present study 
of inert pollutant diffusion and no discussion is given here. However, these 
two submodels are based on very intuitive assumptions coupled with different 
empirical formulations for different physical/chemical processes (e.g., 
descriptions of the aerosol growth and sedimentation velocity are both derived 
from the empirical studies of water droplets in natural clouds). 


3. TREATS Model 

An integrated moment scheme (over the entire horizontal plane at 
each altitude) is used in the TREATS 13 model to determine the dispersion param- 
eters 0^,0 (related to second moments). A Gaussian distribution assumption 
is then adopted to describe the concentration distribution on a horizontal 
plane for each altitude. This model starts with the mean concentration diffu- 
sion equation which is closed by using the diagonal eddy dif fusivities , i.e., 

/ 3 , 3 3 - fiv ixi? Ixlr _JL\ r 

\ 3t + u 3x + V By + W 3z ) C V 3x ** 3x + By By + 3z Kz 3z j C 


+ X C 


( 10 .) 
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where C is the mean concentration of contaminant, u, v, w are the x, y, z 
components of the mean velocity, K x , Ky, K z are the eddy dif fusivities and X 

is a decay constant (which is set to zero in this study) . A four-step proce- 
dure is followed: (1) both sides of Eq. (10) are multiplied by x n y m , where 
0 <n + m< 2 ; ( 2 ) the equation is then integrated over the horizontal plane 
at a fixed height; (3) the boundary assumptions // C(x, y, z, t) dxdy = 0 

X -*■«> 

or y +°° 

and incompressibility assumption 3u/3x + 3 v/ 3 y = 0 are applied to the integral 
equations from (2) ; and (4) the dif fusivities are assumed independent of the 
horizontal coordinates. The following set of equations is obtained; 
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All the quantities in the above equations depend only on the vertical co- 
ordinate z and time t. x(z,t) is the maximum center location of instantaneous 
concentration, i.e., 

x(z,t) = 


where 


nm 


a 


x n y m C(x,y,z ,t)dxdy 


0 oo is the total mass of pollutant at altitude z and time t. 

If the eddy dif fusivities Ky, K z are somehow known, Eqs. (11) 
to (14) form an initial value problem of a set of differential equations 
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which can be solved numerically • A standard finite difference scheme is 
used in the TREATS model. Evaluation of the eddy diffusivities K y and 
K z , which are not known in general, proceeds as follows: 


a. Horizontal Eddy Diffusivity K v and K y - The method incorporated 
to determine the eddy diffusivity is based on the assumption that 
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The standard deviation a of the crosswind concentration distribution is given 
by a functional expression which depends on the local stability. The functional 
expression is derived by fitting the Pasquill-Gifford 29 curves for the corre- 
sponding stability class. The local stability class at each altitude is 
assigned on the basis of the local temperature gradient. The alongwind eddy 
diffusivity K x is assumed to be equivalent to in the present study. 


b. Vertical Eddy Diffusivity K z - Two different methods of 
evaluating K z are used in this study for TREATS. One follows the same pro- 
cedure as that for Ky; the Pasquill-Gifford curves for the vertical disper- 
sion parameters o z are utilized. The other alternative approach adopts 
Blackadar's formulation 30 for the eddy viscosity. If the turbulent Schmidt 
number is assumed to be one, the eddy diffusivity K z will be 


k, ■ *■}(£)* * (g)T / * ♦- 

where <j> is the non-dimensional wind shear and £ is a mixing length. (Blacka- 
dar's formulation is based only on mixing length theory.) The non- 
dimensional wind shear c|> can be derived from similarity theory. In this 
study, expressions for c() have been based on the studies of the surface layer 
for unstable, stable and neutral atmospheres (see Section IV. B). Two kinds 
of mixing length expressions are incorporated. They are 


k(z+z 0 ) 
k(z-fzp) 
1 + A 


( 15 ) 
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where 




2.7 x 10 -4 V 

s. 

f 

if x lcz f 

and £ = 0.0063 — 2- tanh ■ - - ■ (16) 

f 0,0063 u*o 

In these equations V_ is the geostrophic wind speed, k is the Von Karman 

O 

constant, f is the Coriolis parameter, and u and v are, respectively, the east 
and north components of the wind at the level where K z is being evaluated. 

Expression (16) 31 for mixing length £ was obtained by fitting 
the measurement data collected from towers ranging from 18 to 150 m high at 
KSC, where similarity theory was used. 


4. Atmospheric Diffusion, Ihrticl e-in-Cell (ADPIC) Model 

The simplified mean concentration diffusion Eq. (10) can be rewritten 
for an incompressible flow field in the form; 



(17) 


where ' ' denotes the ensemble average, denotes a vector, K is a diagonal 

~ ~ KVC 

3x3 matrix, and U„ = U ;r— is called the pseudo-transport velocity. 

y C 

Equation (17) represents the conservation of mass of a contaminant 
material C in a fictitious flow field of velocity Up. Imposing a strong assump- 
tion that the fictitious flow field is of constant density (or negligible expan- 
sion rates), i.e., V«Up = 0, Eq. (17) implies that a contaminant element will 
travel following the pseudovelocity Up. The ADPIC model 14 is developed on this 
basis. The numerical procedure can be described as follows: 

1. The flow field is divided into a number of cells. Based on the 

initial concentration distribution and the computer storage used, 

each cell is given a number of discrete contaminant elements. 

~ KVC 

2. Using each cell as a unit, the velocity term Up = - — — due to the 

W 

turbulent flux (diffusion), which is the second term of pseudo- 
velocity, Up, is calculated by a finite difference scheme. Up is 
then added to the actual wind velocity U at each cell corner to 
yield a pseudovelocity Up. 
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3 . 


Each tagged contaminant element in a given cell is transported for 
one time step At with velocity Up which is interpolated from the 
pseudovelocity Up at the corners of the cell. In other words, the 
tagged element is relocated in a new position 5^^ given by 


x 


= x 


+ U «At 
P 


new old 

The new concentration distribution after a time step is thus obtain- 
ed on the basis of the new position of each element. 

The eddy dif fusivities 1^, Ky, K z currently incorporated in the 
ADPIC code for this study take the following forms: 


a. Horizontal Dif fusivities - Using the Kolmogoroff theory 32 in 
the inertial subrange of the turbulence eddies and the isotropy ass um ption in 
the horizontal plane, the horizontal eddy diffusivity can be written as 

Kjj = Ky oc e 1 / 3 W 3 (18) 

The proportional constant in the Obukhov expression (Eq . (18)) is taken as one 
and the length scale, £, is given by 

£ = a x (t) = (a 0 2/3 + | £ 1 / 3 t) 3 / 2 (19) 


where a 0 is the initial standard deviation. The rate of energy dissipation, 

£, is assigned as constant (= 2.0) in the model. 

When Eq. (19) for length scale reaches its maximum at long 

times, a slight modification for K , K , based on the work of Walton, 33 is 

x y 

used, i.e. , 


K 

x,y 


(t) 




_1 

K 

max 


+ 


K' (t) 

x,y 


where the K Xj y(t) on the right hand side is that obtained from Eqs. (19) and 

(18) and K is a constant (e.g., K „„ = 5 x 10 9 cm 2 sec -1 for Oo = 204 m) . 
max max 


b. Vertical Ed dy Diffusivity - Since the turbulent shear stress 
in the surface layer is approximately constant, the eddy viscosity can be 
obtained from 

\ < 20 > 
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would be a 


From similarity law, the nondimensional wind shear <j> = — 

function of nondimensional altitude, z/L, where L is the Monin-Obukhov length 
scale (see Section IV.B.l). The functional form for <j> in different stability 
conditions is well known; hence the eddy viscosity can be put into the form 

% - =(t) i (21) 

Using the Reynolds analogy gives Kjjj = K z from the above expression for the 
eddy diffusivity. 

Above the surface layer ( ^ 75 m in the ADPIC model) , K z is 
taken to be constant at its value at the 75 m calculated from Eq. (21). It 
should be noted that the widely used assumption that the vertical eddy diffusi- 
ivity is constant above the surface layer is arbitrary. 


5. A Model for Diffusion in Shear Flow (DISF) 


The DISF model 15 is obtained by analytically solving the diffusion 
equations via a second-order eddy diffusivity closure scheme and the 
Lagrangian approach for statistical quantities. The flow field is assumed 
to be a uniform shear gradient layer. 

Neglecting the molecular diffusion effects allows one to write the 
governing differential equations in Eulerian form for the concentration field 
as 


9C 9C 

3t u± axi 


0 


( 22 ) 


The assumption of indelibility of the tagged points enables identification 
of the instantaneous point source at the origin as the initial condition, 
i. e. , 


C(x,0) = 5(x) 


(23) 


The appropriate solution to Eq. (22) with the initial condition (23) is 

C(x,t) = 6(x - x(0,t)) (24) 

where x(0,t) is the instantaneous position of the tagged element at time t. 
This solution is averaged over the whole ensemble space to give the expres- 
sion for the mean concentration of contaminant , 
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C(x,t) = P(5; 0,t) (25) 

where P(x; 0,t) is the probability density distribution function for x(0,t). 

The dispersion tensor (of which the dispersion parameters 

a x> a z are the diagonal terms) can then be written as 

a ±j (t) = j X i (0,t)X j (0,t) P (x; 0,t) dX C(i,t) x ± (6,t) Xj (6,t) dx 

(26) 


The differential equation for the mean concentration C(x,t) can be derived 
directly from Eq. (22). If the eddy diffusivity tensor is used to close the 
turbulent flux term, the equation becomes 


3 

3t 


+ 



(£,t) 


3x h 


K_ (x,t) 


3 


3x . 
3 


C (x,t) 


(27) 


Equation (27) is the general form of Eq. (10) where only the diagonal terms 
of the eddy diffusivity tensor are considered. 

To obtain an analytical solution of Eq. (27), the eddy diffusivity 
tensor is assumed to be space-independent. Based on Eqs. (25), (26) and (27), 
a set of equations relating the diffusivity to the dispersion tensor 
can be obtained. 3i * That is, for a free flow field with constant shear velocity 


_d_ 

dt 


d 

dt 


2S 

CTi 3 

= 2 K 

- S 

<^3 3 

= k 13 

d 

dt 

O 2 2 

= 2 K; 

d 

dt 

(?3 3 

= 2 K 


with K 23 = K 3a = Ki* = Ka x = 0 (28) 

Furthermore, one can obtain the following expressions by applying 
straightforward Lagrangian methods to the dispersion tensor. 
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VI 2 f Rh(t) dx + S VI VJ f tRi 3 (t) dT 
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Kis + K 3X = VIV^ 


/"[Rax 

J o 


(T) + Rl 3 (x)]dx + S 


J o 


xR 33 (x) dx 


K 22 



(x) dx, K 33 



(x) dx 


( 29 ) 


where R^ is the Lagrangian correlation tensor, is the mixed Eulerian/ 
Lagrangian fluctuation velocity and S is the constant mean shear gradient. 
The eddy diffusivity tensor is non-diagonal as can be seen by examining 
Eq. (29). 

Substituting Eq. (29) into (27), the solution to the mean concen- 
tration equation can be expressed as a generalized Gaussian distribution 


C (x , t ) = — exp I- \ (x - x c ) *A 1 »(x - x 0 )( (30) 

(2xr ) 3 ^ 2 | A | 1 / 2 1 2 ’ 

where A ~ j | » Q is the source strength, x 0 is the source location, and 
superscript T denotes the transposition of a vector. It should be repeated 
that (30) is a solution for an instantaneous point source in a free uniform 
shear flow under an assumption of nonspacially dependent K^j . 

Based on Lagrangian expressions of G^j given by Corrsin, 35 a set 
of explicit and feasible expressions for the dispersion tensor "over all 

the time ranges" is proposed. 


0 2 

= On = 

Cit 2 (Ait + exp (-Bit)) 

X 

V 

— ct 22 = 

C 2 t (1 - exp (-B 2 t) ) 

a z 2 

= g 33 = 

C 3 t (1 - exp (-B 3 t) ) 

a 

X z 

= a !3 = 

Cut 2 (Ae, + exp(-B4t)) 

where the constants are 

given by 
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Ax = | S 2 T 33 t= 4~ 
3 Vi 2 


, Bx = 


Ax - § S 


ViVl 

vi 2 


, Cx = 


Vi 2 


B 2 - 1/ (2T 22 ) 


C 2 — 2T 22 V 2 


B 3 = 1/(2T 3 3 ) 


C3 2T 33 V3 


Aa ST ; 


V£ 2 


B -is Zil 

r>/, — n j „ 


c 4 = VlVl 


ST 33 VI 2 


In these expressions for constants, T 22 and T 33 are the Lagrangian integral time 
scales in vertical and lateral directions. 

The effects of boundaries are accounted for by assuming an 
instantaneous fictitious concentration whose maximum concentration core is 
located at a point x OT The maximum fictitious core Xm is determined through 
a simple analytical geometry calculation by means of an isopleth of 
concentration. 

A relationship between the fictitious core Xoo and the actual core 
x 0 is, for the lower boundary 

(Xoo> yoo, Zoo) 

*0 - \2 z„ ! 

U33 

- Zo 

yco = y 0 


Xoo — 

Xco = 


and for the upper boundary 

Xoo 

Yco = Yo 

Zoo 2H Zo 


= Xo — 


2 (H - Z 0 ) CT x ; 


where H is the thickness of the flow field. When boundaries are appropriately 
accounted for, as indicated above, the instantaneous concentration distribu- 
tion can be written as 
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C(x,t) 


Q 

(2ir) 1 / 2 |Aj 1 / 2 


/ 1 ~ T ~ ^ 

< exp [ — 2~(x - x 0 ) •A - ' 1 • (x - x 0 )] 

+ Y g e xp[-|-(x - Xoo) T -A -1 » (x - Xoo)] 

+ Y u exp[- | (x - 2 Ho ) T .A- 1 (x - ^ o )]} 

(32) 


where y and y are the constants accounting for the absorption strength of 
8 u 

the ground and upper boundaries, respectively. 

The dosage at a position x for a point source, initially at x 0 , can 
then be obtained through integration 


D(x) 


£im -p— / C(x,n)dri 
X \ Jo 


where X is the alongwind distance from x 0 and is the mean wind speed at 
x 0 . Since this is not feasible analytically, particularly with the proposed 
specifications of the dispersion tensor (31) , a Runge Kutta numerical integra- 
tion method is employed. 
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B i INTERPRETATION OF ATMOSPHERIC DATA FOR DIFFUSION STUDY 


The meteorological information required for the diffusion calculation 
is interpreted from the rawinsonde measurements for mean quantities (the meso- 
scale properties) and data gathered at thirteen towers with heights varying 
from 18 to 150 m for turbulent quantities. The towers and instrumentation at 
Cape Kennedy are described in Ref. 36 and will not be described in detail here. 
The meteorological data for each of the Titan III launches have been summar- 
ized in Refs. 37 and 38. The locations and heights of the towers are given in 
Fig. 9. Mean quantities such as wind speed, wind direction, and temperatures 
were available from some towers and at some of the following levels: 1.83, 

3.66, 16.46, 46.34, 62.20, 89.02, 119.82 and 150.91 m. From each tower the 
standard deviation of wind azimuth angle is given at the 1.83 m level; the 
lapse rate, which is the temperature difference between 1.83 and 3.66 m is 
also available. In this section, emphasis is placed on the quantities 
required in calculations using the five models described above. 

1. Atmospheric Turbulence in the FBL 

The planetary boundary layer is a turbulent flow field which gener- 
ally comprises three physical layers: a surface layer near the ground, a mix- 
ing layer at the top, and a free convection layer in between. The surface 
layer is, by definition, that region in which vertical variation of the trans- 
port mechanism characteristics, such as friction velocity and heat flux can 
be ignored. It is usually considered to be 30 m high. However, a recent 
study 39 of the turbulent wind field below 150 m in the KSC area indicates 
that relationships which are valid in the surface layer may apply up to the 
top of the Kennedy tower (150 m) . In addition, a measurement in Minnesota** 0 
also showed that the relationship at the surface layer can be extended to the 
bottom of the mixing layer (approximately 15% of the PBL height). Therefore, 
this study considers the PBL to consist of only two layers, one below and the 
other above 150 m. Interpretation of the turbulence in each layer individu- 
ally is discussed below. 

a. Turbulent Wind Field Below 150 m - The statistics of 
atmospheric flow over homogeneous terrain in equilibrium in the surface 
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layer, and thus up to 150 m are determined by three parameters; the roughness 
length z 0 , the friction velocity, u* and the Monin-Obukhov length, L. 

Based on similarity theory, the diabatic wind profile can be 

written as 

U = — ~ [ln(ze - ^) -ln(z 0 )J (33) 

where k is the Von Karman constant and can be assumed to be 0.4. 

The function ip in Eq. (33) is recommended from previous 
studies 3 9 » ** 1 to be 





2tan 1 ~ — I- tt/2 
Lq 


unstable 

neutral 

stable 


(34) 


where x = (1-16 z/Lo) 1 /** and L 0 is the Monin-Obukhov length at the lower 
portion ( < 30 m) of the layer. 

The u* o and z 0 thus can be determined by using Eqs. (33) and 
(34) if L 0 is known, i.e., u* Q is proportional to the slope of the curve 
plotted from the tower measurement of mean velocity vs. ln(ze~^), and z 0 is 
given by the intercept of the curve with the ordinate. 

The Monin-Obukhov length L is defined as 


L = - 


u * c 
* r 


p T 


k g H 


where is specific heat at constant pressure, T absolute temperature, 
density, g gravity constant and H = p w’q' the vertical heat flux, 
length L 0 is related to the flux Richardson number Rf in the lowest. 10 


(35) 

P 

The 
m by 


L 0 = z/R f (36) 

Since the measurement of w'q ' is not provided, the gradient Richardson 
number is used instead, i.e., 

% - R g - - V ' (ii ) 2 <37) 

where y^ is the lapse rate of temperature and y^ is the adiabatic lapse rate. 


34 





where is the eddy viscosity, is the eddy conductivity and <p is the non- 
dimensional wind shear given by. 


/(l - 16 z/L) -1 / 4 

4> =(l 

^1 + 4.7 z/L 


unstable 

neutral 

stable 


(38) 


The length L 0 is thus determined by applying Eqs. (36) and (37) to observa- 
tions of wind and temperature at 1.83 m and 16.46 m; z in Eq. (36) is given 
as the average value (9.15 m) of these two heights. With the length L 0 , u* q 
and z 0 at each tower have been interpolated using Eqs. (33) and (34) for 
each of the Titan III launches. The observed mean wind velocity used in 
the calculation is based on 30 min averages. In the present study the 
average value of u* o measured at towers along the path of the cloud and at 
the highest towers (110 or 313) is used for each launch case. 

Based on similarity theory, the standard deviation of wind 
fluctuations, a u , a v and a w at alongwind, crosswind and vertical directions, 
respectively, are given as 


a u 

= A u u* 


a v 

= K u * 


CT W 

= Aw u* 

(39) 


The constants A v , are not known precisely. Table VII lists, as a 

reference, the values obtained in different experiments. In this study, the 
values suggested by Yaglom 43 are employed. In fact, the constant expressions 
for A u , Ay, and A w are only true for a near -neutral atmosphere. Experiments 
at Kansas 4 ** and Minnesota** 0 as well as the numerical simulation 45 of an 
unstable PBL suggested that the ratio of <J„/u* can be fit by an expression 

w O 


a 

— — = 1.89 (- z/L) x / 3 (40) 

u*o 
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A more general form has been also proposed in recent literature, i.e. , 

cr 

~ = 1.3(<J> - 2.5 z/L) 1 / 3 (41) 

u * 

where is given in Eq. (38), Although little is known about the ratios 
o u /u A and 0 V / u* o> the ratio O v /u* o is found to be very sensitive to stability. 
Nevertheless, a constant expression for them seems to be best so far. 

In our calculation of 0 W using Eqs. (40) and (41) for the 
August launch in unstable atmospheric conditions below 150 m, the ratio of 
O w /u* below 100 m is within 15% of 1. 07 from Eq.(40) and 1.37 from Eq. (41). 

The average value of these two constants is surprisingly close to Yaglom’s 
value for A , 1.25. 

t.t 7 


b. Turbulent Wind Field Above 150 m - In the limited literature 
dealing with turbulence above 150 m, no consistent results are to be found. 
In this study the work is based on a simple governing equation of a steady- 
state, neutral barotropic PBL, i.e.. 



_d_ 

dz 


u’w' 


£(V - V g ) 


(42) 


where V and V are the component of above ground wind and geostrophic wind, 
respectively, at right angles to the surface wind. 

The geostrophic wind is taken as the wind at the top of the 
PBL. The calculation for u* 2 using Eq. (93) is started at 10% of the height 

of the PBL with the initial value u* q 2 . Figure 10 shows the relative value 

u* 2 /^ 2 vertically for the December and August cases. The constant value 

for \i* in the lower portion is found to have up to 5% error up to 1/6 of the 

PBL. 


In view of the lack of information regarding a u , a v , and o w 
it is assumed that the relationship of Eq. (39) is also valid for the layer 
above 150 m. This seems reasonable for a strongly convective PBL in which 
the mixing layer is nearly neutral. (The results from the Minnesota experi- 
ment support this point.) However, it should be kept in mind that neither 
of the three selected launch cases is in a strongly convective condition. 

A determination of these basic quantities should be obtained if the 
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capability to make confident environmental assessments is desired. 

Since the mean wind profile for the August case (Fig. 3 ) shows 
a negative wind shear between 500 and 1000 m, the positive vertical momentum 
flux in this region would imply that the eddy viscosity is negative, and that 
there is thus a negative eddy diffusivity by the Reynolds analogy. Although 
negative eddy dif fusivities may (and here do) occur occasionally in descrip- 
tions of strongly bounded flows, none of the codes employed can handle them. 
Therefore, a negative value of u* 2 in the region between 500 and 1000 m as shown 
by the dashed line in Fig. 10, is used in the DISF calculation because DISF was 
the only model requiring u* as input . 

2 . Mesoscale Meteorological Properties 

The required mesoscale meteorological information such as wind speed, 
wind direction, temperature, pressure, and relative humidity in the PBL is 
based on the sounding measurements performed by the Air Force Eastern Test 
Range Weather Group. There are three kinds of measurements (rawinsonde, wind- 
sonde and jimsphere) at different times. Since only one measurement is given 
at a time at KSC, a full-scale description for the atmospheric flow at the 
KSC area depends on the theoretical model. In this study, all models except 
ADPIC are based on the assumption of homogeneity in the horizontal plane and 
stationarity during the period of calculations; a simple meteorological model 
to simulate the full scale wind field accounting for the spatial and time 
changes is discussed in Section IV. C. 1 . The present study uses mainly rawin- 
sonde measurements; windsonde and jimsphere measurements are used only for 
reference. The wind speed, wind direction, humidity, and virtual potential 
temperature from rawinsonde measurements for the three Titan III launches , 
(the set of measurements at nearest time of launch is used) , are shown in 
Figs. 2, 3, and 11. 

The virtual potential temperature is calculated using Tabata's 
expression, " 7 i.e.. 


v 


T r 1 + 1 - 61 w m j r iooo i 0,288 
L l + w JLpJ 

m 


where T is air temperature (°K), P barometric pressure (mb) and W m is mixing 
ratio given by 
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0. 622»RH»e 

y = — 

m P - RH«e s 

where EH is the relative humidity (%), e g is saturation vapor pressure = 
(c-dx-ex 2 ) 

10 , x = 1000/T, and the constants c, d, and e are given by 

c = 8.42926604, d = 1.82717843, e = 0.07208271. The virtual potential tempera- 
ture is plotted instead of the air temperature because it is the main param- 
eter for determining the thickness of the PBL, the existence of the inversion 
layer, the stability of the layer, etc. 

As mentioned before, a sensitive and important meteorological 
parameter in the prediction of the cloud ground concentration is the thickness 
of the PBL. The PBL in this study is defined as the region in which the 
turbulence cannot be ignored in the transport process. In a well-mixed PBL 
capped by a temperature inversion, e.g., in sunny (clear) daytime, it will be 
approximately neutrally stratified at the upper portion and its mean wind 
profile will be nearly uniform vertically. Idealized characteristics of the 
well-mixed PBL are shown in Fig. 12 . From the plots in Fig. 12 , it is evi- 
dent that none of the selected launch cases is similar to this ideal condi- 
tion; in fact, the August case is a sea-breeze condition. At Florida it is 
shown that the sea breeze usually initiates a storm;'* 8 thunder during the 
launch time was indeed heard at the KSC weather station. The March case was a 
late afternoon launch (at 2027 EDT) and the December case was a night launch 
(at 0310 EDT) . Since the detailed data on turbulence needed to establish 
the thickness of the PBL is not available in these cases and since an ideal 
PBL did not exist during these launches , the PBL thickness must be inferred 
from observations of cloud behavior, as discussed below. 


a. August Case - As shown in Fig. 3 , the virtual potential 
temperature increases sharply after 1600 m, the wind speed changes dramati- 
cally after 1850 m, and the relative humidity drops significantly between 
400 and 700 m. These observations are not consistent ; any determination based 
on these mean quantities would not be very meaningful. Alternately, we can 
use the simplest balance Eq. (42) and integrate it, and obtain 


u 


2 

* 


(z) 



fJ'J (V g - V)dz 


(43) 
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If the top of the PBL is assumed to be a zero vertical momentum flux, 
u* 2 (Z^) = 0, where Z-^ is the thickness of the PBL, and Vg = — A u* 0 , which is 
the Kazanski-Monin relationship, Eq. (43) becomes 

/* Z i 

u* = — Af Z-u* - / f-V dz (44) 

J o 

Since u* o 2 is known, f is the Coriolis parameter taken to be 
5,5 x 10 -5 (sec -1 ), and V is known at each height, we can obtain Z^ by solving 
Eq. (44) for a given A. It is found that only one Z^ exists for the value of 
A between 5 and 15.^: Z^ is found to be 1800 m. We have found that u* 2 equals 

approximately 0 at 1950 m in Fig. 10, where a different approach was used 
(the geostrophic wind was assumed to be the wind at 2000 m) . The close agree- 
ment between these two results is not surprising because they are both derived 
from the same equation and use very similar values for Vg. The prediction of 
2000 m for the thickness of the PBL does not agree with the value of 1090 m 
given in Ref. 23. It is not clear how that value was derived. However, real- 
time tracking of the cloud by camera shows that the cloud dispersed faster 
below 2000 m and that an inversion cap existed at about 2000 m. In other 
words, observation implies that turbulence cannot be ignored up to 2000 m and 
thus the thickness of the PBL will be no less than that height. 

b. December Case - Using the approach of the August case for the 
December case, the thickness of the PBL is found to be about 510 m. It is 
slightly different at the 540 m height, where the u A 2 is zero, shown in Fig. 

10 . The virtual temperature increases strongly after 620 m and the humidity 
drops significantly at about the same height as shown in Fig. 2, so again the 
thickness of the PBL seems to be about 600 m. In fact, the results from air- 
borne measurements 1 1 at an altitude between 560 m and 610 m demonstrate a very 
slow dilution of the cloud in this region indicating very low turbulence levels. 
Since the Titan ground cloud is always stabilized at pretty high levels (with an 
order of magnitude of 1000 m) , most of the exhaust mass will stay above the PBL 
under a shallow PBL condition. If the PBL is assumed to be thicker, more of the 


t A = 12 is often used in conditions without data. 
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exhaust mass will remain in the PBL. Since increase of mass in the PBL will 
certainly increase the ground level prediction, for conservative purposes we 
will use the higher value 660 m (the height at which measurements were obtained) 
for the PBL thickness in the diffusion calculations for the December case. 

c. March Case — In this late afternoon case, all the meteorological 
parameters have a turning point between 700 and 1100 m (see Fig. 11). The IR 
observation of the cloud showed a separation of the atmospheric flow at about 
1100 m. Since at this time of day one expects to be in a period of transition, 
any estimate will be subject to significant error. For the calculations we 
simply take the thickness as the sounding measurement height of 1088 m. 

It would be interesting to know the actual ratio of the selected 
thickness and the widely used length scale u*/f for the Elcman layer; the 
ratios Z^f/u* used for the August, December, and March cases are 0.18, 0.16, 
0.14, respectively. They are all lower than the value 0.28 obtained from a 
simulation of the idealized neutral atmospheric boundary layer. 49 However, be- 
fore a more precise micrometeorological model is established to provide the 
information on the PBL height, we may use 0.16 u A /f to be a reference height for 
determining the PBL height from the sounding measurements. 

C. CONCENTRATION CALCULATIONS 

Although each model selected for the diffusion calculations uses 
different forms of input data, the basic information required by each is 
(i) the initial cloud shape and concentration distribution, (ii) the advec- 
tive wind field and boundary conditions, (iii) the statistical/physical 
quantities needed (generally different for each model) to describe the turbu- 
lent mechanism, and (iv) the type and location of output desired from the 
calculations , 

This basic information for the aforementioned three Titan III 
launches is given in Appendix C. The cloud shape and concentration distribu- 
tion at cloud stabilization are calculated using the NASA/MDM preprocessor 
program 1 for the layered model 4; the meteorological data and the thickness 
of the PBL are determined as described in the previous section; surface 
measurement data are taken directly from the tower measurements. All the 
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statistical parameters required for calculations by TREATS, ADPIC and DISF 
are interpolated by following the procedure given before. Calculation 
methods using the MDM program have been covered in the parametric studies of 
Section III; diffusion calculations using the other four models are described 
below. 

1. The METS Model 

Utilizing the input information given in Appendix C, calculations 
using METS were made employing several options in the program to provide an 
upper bound of its predictions. The METS calculations that include parti- 
cles and chemistry give a higher ground prediction of HC1 than if they are 
not included. The particle size distribution in the calculation is chosen 
randomly . 


2. The TREATS Model 

In this study, the flow field is divided into 45 sublayers, at each 
of which moment integrations over the entire horizontal space are calculated. 

The source strength and input information for each sublayer are obtained by 
linearly interpolating values provided in Appendix C. The time step for the 
calculation is set at 30 seconds. 

For the August study, three calculations have been made using differ- 
ent eddy dif fusivities K z as given in the review section of the TREATS code. 
Since the value of K calculated from the Blackadar formulation (15) and (16) 
is about one and is probably too small for z ^ 100 m in a convective PBL in 
general, a more probable value of K z = 30 for z ^ 100 m is used in the August 
calculations. The value 30 is calculated using the expression for K z based on 
the similarity theory as given in the review section on ADPIC. There is only 
one calculation for the March study and two for the December study. The March 
calculation uses the K z from Pasquill-Gif ford curves and the December calcula- 
tions use the K z from Pasquill-Gif ford curves and the Blackadar formulation (15) 
and (16). The parameters required for the calculation in the Blackadar formu- 
lation, (e.g., z 0 and u* q ) use, for all three launches, the average value of 
those interpolated from tower data using similarity theory. 
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3. 


The ADPIC Model 


Because of limits imposed by the availability of computer storage, the 
total number of grids for each calculation is assigned to be 41 x 41 x 15 = 
25215. The grid size is 150 m vertically, 500 m laterally, and 500 m longitu- 
dinally in the August study; 80, 500, and 500, respectively, in the March study; 
and 50, 1000, and 1000, respectively, in the December study. The maximum number 
of source terms is five. Each may be located at an (x, y, z) with a Gaussian 
distribution of material in three dimensions. Each dimension has a right and 
left (or up and down) cutoff distance- from x,y,z. Each source may also be a 
different species with individual particle distribution. The total number of 
species plus locations must be less than five. Initially, the concentration in 
each pancake is assumed to be uniform vertically and Gaussian in the horizontal. 
Each sample point is assigned a weight of 819 g, 444 g and 135 g for August, 
March and December, respectively. The total source is represented by 20,000 
samples in each case. Each cell has many more than one sample particle. 

As described in the brief review of the ADPIC program in the 
previous section, the program can utilize any time averaged advective wind 
field as input at a given time. Therefore, the ADPIC program can account for 
the variation of wind field downstream. Since only one mean wind measurement 
(from rawinsonde) at KSC is given at intervals of two hours or longer, the 
change of wind field with respect to space as well as time must be estimated 
if the nonhomogeneity or nonstationarity of the atmospheric wind field is 
to be considered. A simple three dimensional wind field model 50 is employed; 
it utilizes the surface tower measurements because they provide the space and 
time variations. The procedure to create a wind field is as follows: (1) It 

is assumed that the vertical wind at each tower location has the same profile 
as the rawinsonde measurement but different wind direction and magnitude of 
wind speed; the wind direction and the magnitude of the wind speed are deter- 
mined by fitting the lower portion of the assigned profile to the measured 
wind data of the tower. (2) The wind field at each grid point is then inter- 
polated from the wind value at the tower locations by allowing a minimal 
adjustment which relies on the constraint of incompressibility of the flow 
field (i.e., V»u = 0). (3) The simulated wind field is changed whenever 

the tower measurements change. Each of the sixteen towers shown in Fig. 9 
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is used in estimating the wind field. Two sets of wind data, separated by 
30 minutes each, are used in the August study; ten sets and twelve sets over 
five minutes each are used in the March and December studies, respectively. 


4. T he DISF Model 

Since this model is formulated for a point source, a diffusion 
calculation for an irregular non-point source would require an integration 
of point source calculations over the entire initial cloud, i.e.. 


C(x,t) 



x 0 , 


t) dx 0 


(45) 


where x 0 is a point position in the cloud and C(xJ x 0 , t) is the concen- 
tration at x at time t from a source point Xo, as given in Eq. (32). This 
is not practical because it cannot be integrated analytically and numerical 
integration requires large computer times. An alternative approach is to 
divide the initial cloud discretely into a number of cell sources, each of 
which is treated as a point source. The total number of cell sources is 
assigned to be 40 for each layered pancake in the present study. The number 
of pancakes is the same as that given in Appendix C. In order to reduce the 
error induced by the discrete approach, each layer is divided in the follow- 
ing way, which shows a cross section 



where O is the standard deviation of concentrations at the time of cloud 
stabilization. The point source is located at the center of each cell. The 
physical quantities used (e.g., mean wind field) for a calculation from a 
point source are determined by the following method: 
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1. Wind direction is the average value of wind directions in the 
range of altitude between 1/6 of the PBL thickness downward from 
the altitude of the point source and 1/6 of the PBL thickness 
upward. Wind directions outside the boundaries are ignored. 

2. The wind shear and wind speed are derived from the best linear fit 
of the measured wind data in the range cited above. 

3. The turbulence intensities and stress are assumed to be constant 
below the 1/6 depth of the PBL at the value interpolated from tower 
data using similarity theory (see Section IV. B). Above 1/6 of the 
PBL, the average value of an assigned turbulence profile within 
the range given in (1) is used. 

D, COMPARISON 

1. Model to Model Comparison 

The results calculated from all five models for the three Titan 
launches are shown in Figs. 13 to 19. Three kinds of values from the calcula- 
tions are plotted for the comparison. They are (1) maximum instantaneous con- 
centration, (2) maximum integrated concentration (dosage), and (3) the path of 
maximum ground level concentration. Only one species (HC1) is considered. As 
mentioned before, some of the models used to make calculations employ estimated 
values for some input parameters which may not have much physical validity. At 
present, only the plot that gives the highest value for each model (except MDM) 
is given in each comparison. There are at least two plots from MDM calculations 
in each figure showing concentration or dosage. The reasons for doing this are: 
(1) Although it was generally recommended that the MDM calculations should be 
performed using tower measurements of O A , we used both tower data and the inter- 
polated 0 Aq from similarity theory (which are more consistent with the input 
requirements of the other models). (2) It is of interest to see whether the un- 
certainties from a basic input parameter, 0 Aq , which is used to describe the 

turbulence mechanism in the MDM will cover the range of predicted results from 
other models. (3) We hoped to determine the best approach for using MDM (in 

terms of doing an environmental impact study or determining launch constraints 
should the MDM be used for these purposes in the future) . 
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Since the use of the TREATS model with the Blackadar mixing length 
formulation (15) and (16) yields the highest value, the plots in Figs. 13 to 17 
represent results using this approach. The dosage plots in Figs. 14 and 17 
from the METS model were obtained from calculations which include the chemistry 
and particles because otherwise near zero values would result. 

The figures and comparisons are discussed below for each launch case. 

a. 20 August 1975 Case - For this typical afternoon sea breeze 
weather condition at KSC, ADPIC, TREATS and DISF predict nearly equal maximum 
instantaneous ground concentrations downstream, especially in the near field of 
the launch complex, while results from the MDM give a higher level. Specifi- 
cally, the plotted results in Figs. 13 and 14, from the MDM model 3 and model 4 
calculation using the value of cr^ given in Ref. 23, show a very high maximum 
concentration in the near field although in the far field they show a better 
value, nearly equal to those obtained from the other three models. It should 
be noted that all the predictions will converge far downstream ( *» 100 km) if 
the wind field is frozen. There is a slight difference between dosage compari- 
sons and the instantaneous concentration comparisons. In Fig. 14, it is 
clearly seen that DISF gives higher predictions of maximum dosage than ADPIC 
and TREATS but merges with the MDM results about 30 km from the launch pad. The 
METS model did not provide enough information for a clear comparison, but appar- 
ently predicts values too low in the near field. For a reference, the results 
from TREATS Using two alternatives for the vertical eddy diffusivity (see pre- 
vious section) are in the bound of the 30% variation for maximum concentration 
and 50% for dosage from the plotted values in the figures. It should be men- 
tioned that the plots of Figs. 13 and 14 ignore the fact that each model pre- 
dicts a different transport direction at downwind distances. Figure 15 thus 
gives ground paths of the maximum concentration downstream from the five models. 
In this figure, DISF and MDM have the same path 8 km from the launch complex. 

The TREATS path is the furthest from the MDM and the ADPIC and DISF paths lie 
between the TREATS and MDM paths. It is interesting that the METS model pre- 
dicts a path jumping between the paths of MDM and TREATS. 
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b. 10 December 1974 Case - For this nocturnal case, the 

comparisons. Figs. 16 and 17, show characteristics similar to the 20 August 
1975 runs, except for the following points: (1) All the predictions for the 

maximum instantaneous concentration converge much faster downstream. (2) There 
is a big difference in the near field between the MDM predicted results using 
measured and interpolated a. ; the farmer is much lower than the latter. (3) ADPIC 
gives very low predictions. (4) ADPIC, DISF, and TREATS predict very differ- 
ent maximum concentrations in the near field. (5) The METS model gives a 

low dosage prediction, especially downstream. 

In Fig. 16 for the comparison of maximum concentrations, two 
additional results are shown which were calculated using MDM models 3 and 4 
with an input provided by a different version of the NASA preprocessor pro- 
gram. Clearly they produce very high predictions which are probably unreal- 
istic (see Section IV. F) . 

One additional plotted result from MDM is also given in Fig. 17 
for the comparison of maximum dosage predictions. This line results from the 
use of a vertically varied cf^ for input, as mentioned in Section III.A.l. 

This predicted result seems to correlate better with the results calculated 
from the other models. However, its lower prediction in the near field raises 

some concern that the MDM model may underestimate the ground level concentra- 
tion in the near field ( ^ 10 km) from a rocket launch if better and more accu- 
rate turbulence information for the PEL is available and used in the MDM 
calculation. 

c. 14 March 1976 Case - This is a late evening case. Figure 18 
shows good agreement for maximum concentration predictions from all models 

in the near field except the MDM using interpolated Oa 0 which gives a slightly 
higher prediction on the ground. The dosage plot. Fig. 19, demonstrates the 
much higher prediction from DISF than from all the others after 15 km. This 

high predicted value from DISF downstream is caused by the assumption of a 
vertically exponential decay of turbulence intensity above the surface layer, 
while actually the lower atmosphere at KSC at launch time somehow still main- 
tained some turbulence strength. The existence of strong turbulence above 
the surface layer is shown by the observation that the south cloud (with which 
we are concerned here) dissipated 15 minutes (6 km downwind) after launch. 


46 



Hence it is reasonable to assume that, after moving 10 km downwind, the cloud 
is vertically well-mixed and further dilution is controlled by horizontal 
dispersion. Since the horizontal dispersion rate is increased after 10 
minutes, the maximum dosage should decrease downwind. Therefore, the 
increased profile of the maximum dosage predicted by DISF must result from 
an insufficient vertical dilution at earlier times (or a slow diffusion) . 

Since turbulence is the only mechanism that controls the dilution it is clear 
that the turbulence input for the DISF model is too weak. However, the impor- 
tant advantage of the plotted results from DISF for the present study is that 
they provide a conservative prediction of possible ground level concentration 
far downstream. Since the DISF predicted values are still lower than the peak 
value predicted by MDM using interpolated O. (although this peak value is at 
the near field of the launch pad) the MDM predicted peak of maximum ground 
level concentration seems to be a conservative input for analysis of the ground 
level environmental hazard. 

It should be mentioned that all the results calculated from DISF 
using any profile other than exponential decay for turbulence above the surface 
layer fall in the range below the conservative prediction line drawn in Fig. 19 

2 . Model to Measurements Comparison 

In order to better evaluate the models, their predicted results and 
available measurements have been compared. These comparisons utilize both 
ground-based and in-cloud airborne measurements and emphasize the August and 
December launches. The flight path of the aircraft in the August study, as 
shown in Figs. 20 and 21 is plotted from ground-based radar tracking data but 
in the December study it is based on the airplane crew's visual record of the 
aircraft location during measurements. The HC1 data were obtained from the 
chemiluminescence monitors; detailed information on the monitors can be found 
in Ref. 51. 

a. August 20 Case - Airborne Measurements - Fifteen aircraft passes 
were made through the cloud at altitudes from about 1100 m to 1600 m (see 
Fig. 21). For the sake of simplicity we have chosen six data sets for the 
comparative analysis, viz., passes 5, 6, 9, 10, 12, and 14. The altitudes of 
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these passes for the model calculations are 1420 m, 1125 m, and 1600 m for 
passes 5 and 6, passes 9 and 10, and passes 12 and 14, respectively. 

Since neither the MDM or METS programs currently provide 
information on instantaneous concentrations , direct comparisons between the 
predicted values and the airborne measurements were made only for the ADPIC, 
TREATS and DISF models. Figures 22 and 23 show the comparative results. The 
solid lines in the figures are the values calculated from the model on the 
path as shown in Figs. 20 and 21 at the appropriate time. The dotted lines 
are the experimental results. Because many uncertainties are involved in 
both the calculations and the measurements, it is difficult to derive any 
conclusions as to the validity of the absolute values at any spatial location 
in either case; for example, the wind velocity used in the model calculation 
is certainly not the real wind velocity which transports the ground cloud at 
the launch. This uncertainty could cause the model-predicted cloud to be far 
from the real cloud, i.e., the path of measurements is not the same as that 
simulated. Therefore two additional calculations for each pass have been made 
by simply moving the whole path of each pass 500 m downwind and 500 m up- 
wind. The results are drawn in dashed lines for the upwind relocation and in 
dashed-dotted lines for the downwind case in Figs. 22 and 23, The data bars 
on the solid lines from TREATS represent the range of variation obtained 
using different vertical eddy dif fusivities as given in the previous section. 
Surprisingly, the values fall in a quite narrow band. 

Comparisons of concentrations on flight passes between the air- 
borne measurement and TREATS model are given in Fig. 22(a) and comparisons for 
the ADPIC model are given in Fig. 22(b). The comparisons for TREATS on 
passes 9 and 10 are not plotted in the figure because this model produces a 
value of 0 ppm on all three paths for both passes, i.e., the "exact" 
path, the path moved 500 m downstream, and the path moved 500 m upstream. 

This is not surprising in view of the comparisons for other passes. From 
the plots for passes 5, 6, 12, and 14 in Fig. 22(a), it may be seen that the 
maximum concentration in the TREATS-simulated clouds during the time of each 
flight pass is much higher than the measurement. As mentioned before, these 
four passes are presumed to have traversed the maximum concentration region 
of the cloud. Therefore, the higher values obtained from the model imply 
that diffusion in the TREATS model is not as fast as in the real cloud. In 
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other words, the eddy dif fusivity used in the TREATS model may be inaccurate. 

The slower rate of diffusion prevents the cloud from becoming large enough 
to produce any concentration for passes far away from the cloud center such 
as passes 9 and 10. The comparisons for ADPIC results shown in Fig. 22(b), 
imply the same conclusion as for the TREATS model. No exact calculations 
were made for passes 5, 6 , 9, and 10. 

As shown in Fig. 23 > concentrations calculated from DISF 
correlate, in general, quite well with flight measurements. On the basis of 
comparisons between calculations and data for passes 5 and 6, we can conclude 
that the cloud is actually located a little further downwind than the model 
indicates. This apparent dislocation, however, raises the question of what 
the maximum predicted concentrations from DISF would be on a path through or 
near the center of the modeled cloud where both passes 5 and 6 are presumed 
to lie. The answer is that the maximum predicted concentration anywhere in 
the cloud will be very close to (although definitely somewhat higher than) 
the value of the dashed lines in Fig. 23. This conclusion is based on the 
observation that although the exact location of any of the calculated curves 
with respect to the cloud coordinates is unknown, that position can be 
reasonably estimated by assuming a bell-like spatial distribution at any point 
in time and then comparing the relative concentrations of the three profiles 
which are separated by 500 m. Thus, in pass 5, it may be concluded that the 
"exact" profile (solid line) and the upwind profile (dashed line) both lie 
fairly close to the center of the cloud since this is the only portion of the 
cloud where fairly high concentrations can be encountered and a 500 m shift 
in position does not change the maximum value significantly. Similarly, in 
pass 6 three curves appear to be on the upwind side of the centerline profile 
but, since the differences between them are small, the two highest curves, 
at least, are not far from the center. The 3 ppm maximum in the profile far- 
thest upwind therefore is a fair approximation of the absolute maximum 
attainable. 

On the other hand, passes 9 and 10 exhibit very low concentra- 
tion profiles and display a monotonic increase in maximum concentration as 
cuts are taken further and further upwind. This indicates that this region 
of the cloud lies near its edge where the concentration vs. time distribution 
is characterized by the flat tail portion of the curves and low absolute values. 

49 



Considering the comparisons of flight passes 12 and 14 in 
Fig. 23 , it is found that the peak value of downwind (dash-dot) profile is 
higher than that for the solid and dashed lines, directly opposite to the 
results of the other passes in the figure. This is because the wind speed 
used in the DISF calculation at the altitude of these two flight passes is 
the average value for the upper 700 m of the PBL and the actual wind speed at 
1600 m is higher than the average. The argument concerning the differences 
between plotted results for passes 5, 6, 9, and 10 is not applicable to 
passes 12 and 14 because the separation distance between the two paths in 
the calculations (500 m) is smaller than the standard deviation of the cloud 
spreading when the cloud is better mixed after long times; the better mixing 
(i.e., wider spread) is apparent from the data (dotted lines) for passes 12 
and 14. 

There is an additional interesting point in the ADPIC compari- 
sons for pass 12 in Fig. 22(b), this plot gives an excellent example sup- 
porting the argument used above in the DISF flight comparisons to qualita- 
tively determine the maximum concentration via analysis of the differences 
of the peak values of the plotted results. As shown in the plot, the 
calculated maxima on the given three paths increase rapidly and progressively, 
in the upwind direction. Using the argument cited before, the maximum value 
at the centerline of the cloud should be located even further upwind and its 
value will be significantly higher than the maximum in the dashed curve — how 
much cannot be accurately determined. In fact, detailed calculations of the 
entire cloud using ADPIC yield a centerline maximum value of 10 ppm. 

b. August 20 Case - Ground-Based Measurements - For the August 
launch there were five ground-based HC1 monitoring sites indicated by black 
spots in Fig. 15 . No HC1 was detected at any of these sites except P-10 near 
the coast line. "No HC1" means that the HC1 level is less than the lower 
detection limit, 0.005 ppm, of the instrument. As shown in Fig. 15 , the 
ground level calculations using TREATS and ADPIC predicted that the path of 
maximum concentration would be closer to P-2 , P-5 , P-7 , and P-8 where no HC1 
was detected than to site P-10; TREATS, in particular, predicts a path that is 
farthest from P-10. If a direct comparison is made between the measurements and 
the results calculated for these sites by ADPIC and TREATS, great disagreement 
is seen; ADPIC did not give a larger value than the instrument limit on sites 
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P-7, P-8, and P-10 but gave measurable values at sites P-2 and P-5; TREATS 

gave no HC1 on site P-10 but gave measurable values at the other four sites. 

In order to test these models with less demanding requirements, 

one can ignore the direction of cloud travel and simply assume that all paths 
of maximum concentration for different models coincide along the path predicted 
by MDM because the MDM (as well as DISF) predicted path is closest to P-10. The 
comparison of results for site P-10 under such circumstances is shown in Table 
VIII. DISF gives a value near that of the observed maximum concentration but a 
value which is a factor of 5.6 higher than the measurement for the dosage. 

In contrast to DISF, the TREATS model showed an excellent agree- 
ment with the observed dosage value but a factor of two lower than the measured 
maximum concentrations. Compared to the data, ADPIC predicted too low and MDM 
predicted conservatively high on both concentration and dosage. It should also 
be mentioned that the MDM model is the only one of the four models (excluding 
METS) which predicted the HC1 value larger than 0.005 ppm at those "no HC1" 
detected sites with these adjustments; in fact, it is quite high (about 0.2 ppm). 

c. December 10 Case - Airborne Measurement - During this launch, 
there were 12 sampling flight passes through the lower cloud below 660 m, the 
PLB height in this case. Most were at altitudes between 550 m and 600 m. From 
the flight measurements, the portion of the cloud between altitudes 550 m and 
600 m was located between 180° to 185° from the launch pad. The actual path 
is indicated by the shaded area in Fig. 24 and is about 10° to 15° off the 
cloud paths predicted by the models as indicated by the dashed and solid lines 
in the same figure. The predicted cloud path was calculated using wind direc- 
tion measured at 40 min before launch (see Appendix C) . The lack of correla- 
tion between the measurements and predictions may be attributed solely to the 
change in the wind. Furthermore, the cloud path at the 1400 m level, observ- 
ed using IR measurements and the dashed line in Fig. 24 , is also found to have 
a 10° to 15° difference from the wind direction used in the calculation at 
that 1400 m level. Therefore, we have compared the predicted concentration 
in the maximum core region directly to the measured airborne sampling data. 

The maximum core region was taken as the disc-like volume between altitudes 
550 m and 600 m and with radius 550 m. Its volume is approximately equivalent 
to the cell size used in the ADPIC model. To make such comparisons, we have 
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assumed that the airplane average speed was 55 m sec -1 and the concentration of 
the core region was well mixed, i.e. , the average value of the concentrations 
within the 22 sec interval during which the highest levels were measured was 
used to represent the measurement value of the core region at the time of the 
flight pass. Surprisingly, it was found that the measured value at 28.52 min 
after launch was a factor of three greater than the initial value used in the 
model calculations at stabilization time. This meant that either the total 
source strength assigned to the initial cloud (below the PBL height) using the 
preprocessor program was less than 30% of the strength of the real cloud, or 
the lower portion of the cloud near the ground was unrealistically given too 
much strength. The latter situation is impossible because the total source 
strength of the portion of the cloud in the calculation below 600 m was less 
than three times the assigned source strength of the portion of the cloud be- 
tween 550 m and 600 m. In other words, the preprocessor program underestimates 
the total amount of pollutants left in the nocturnal PBL. This probably occurs 
whenever the predicted cloud stabilization height is above the estimated PBL. 

Further comparisons between the predictions and measurement 
data were made for the so-called "dilution ratio" that is, the ratio of the 
average concentration in the core region at time t to the concentration at 
time t-At. Table IX shows the ratios calculated from different models and 
airborne measurements. The earlier time refers to the times at stabiliza- 
tion and 10 min after. The later time refers to 26 min and 36 min after 
stabilization; 26 min after stabilization is the time of the first airborne 
sample pass made between altitudes 550 m and 600 m. This table shows that: 

(1) ADPIC and DISF models dilute the cloud very slowly at the upper portion 
of the PBL, that is, they move only small amounts of pollutants from the 
upper highly concentrated region to the ground; (2) TREATS seems to predict 
slightly faster diffusion than the real situation; and (3) MDM definitely 
diffuses too fast during the first 10 min at the upper level. The fast diffu- 
sion predicted by MDM at the earlier time explains why the highest ground 
concentration resulted at the near field of the launch pad using MDM. Most 
importantly, it shows that the high peak of maximum ground concentrations 
predicted by MDM is caused by moving the upper level pollutants to the ground 
too fast. It should be noted that the results shown in Table IX from MDM 
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were based on calculations using the vertically uniform a ^ and a E<> (see 
Section III.A.l). 

d. December 10 Case - Ground-Based Measurements - Locations of 
the four ground monitoring sites and the predicted paths of ground maximum 
concentration are given in Fig. 25 . The measurements have shown that no HC1 
was detected at site P-1 but 19.5 (ppm-sec) , 6.2 and 15.2 were collected at 
sites P-2, P-3, and P-4, respectively. With this information, the ADPIC pro- 
gram clearly demonstrates its capability to predict the maximum ground concen- 
tration path of pollutants more accurately than the other models. This is 
probably due to its ability to model the advection effect; the December 
case was in a nocturnal condition with strong advective wind and weak turbu- 
lence. In the comparisons made for the December case, we assumed that all 
paths of maximum ground concentration for different models coincided in one 
path given by ADPIC and then compared the concentrations at locations P-2, 

P-3, and P-4 between the model predictions and the measurements. The 
results are shown in Table X.. Since in the airborne measurement comparisons, 
the thtal source strength of the cloud used in the calculations is about a 
factor of three less than that of the real cloud, the values in the table 
have been multiplied by three. 

We have also given three sets of values for MDM in the compari- 
sons. They were calculated using the vertically uniform = 8°, = 4° and a 

monotonic decay profile for a ^ as cited in Section III, respectively. Appar- 
ently, none of these three sets of values gave trends similar to those observed. 
This may be attributed to the improper dispersion by MDM at the upper portion 
of the PBL. The ADPIC results are clearly lower than the measured data. The 
paradox of why ADPIC gives lower predictions of ground instantaneous concentra- 
tion than DISF, while the former has a greater dispersion than the latter in 
the upper portion (as shown in the airborne comparison) , can be resolved in two 
ways: (1) DISF assumes higher turbulence than ADPIC in the middle portion of 

the PBL and (2) the imposed assumption of incompressibility of the fictitious 
flow field, which would reduce the diffusion rate, in the ADPIC model may be 
improper. Although the maximum instantaneous concentrations predicted by 
TREATS and DISF were notably lower than the measurements at sites P-2 and P-4, 
the predictions from these two models seem to be adequate. However, from the 


53 


overall view of the prediction capabilities of these five models by comparison 
with data, TREATS gives the best results for this nocturnal launch case. This 
may be because (1) since the shear of the wind speed in the PBL is negligible 
in this night case, ignoring off-diagonal eddy diffusivity terms in TREATS is 
valid, (2) the Blackadar type vertical eddy diffusivity and the mixing length 
scale are proper for the shallow PBL as in the December case, and (3) the empir- 
ical Pasquill-Gifford type curves were obtained from measurements at the lower 
portion of the layers; when the layer is shallow, the eddy diffusivity derived 
from these curves may be adequate. 

E, CALCULATIONS AND COMPARISONS FOR AN OBSERVED GROUND CLOUD 

The comparative studies described in the previous sections were 
made utilizing as input the description of the stabilized cloud provided by 
the MDM preprocessor program. (For simplification, we will call it the 
"modeled cloud".) As mentioned before, the uncertainties included in such 
input may seriously affect the downwind predictions of concentration distri- 
butions. Hence a comparative study of diffusion models using more realistic 
input was desired. Of the three launch cases examined in this study, the 
August case is the only one for which the most complete observation and air- 
borne sampling data are available. The most detailed comparison based on 
stabilized cloud characteristics as interpreted from both observation data 
and airborne measurements (collectively called the observed cloud data) was 
therefore made for the August case. Because the previous comparative results 
for this case showed that the DISF model yielded the best predictions com- 
pared to measurements, this set of calculations and comparisons for the 
observed cloud focused on the DISF and MDM models only. The assumptions and 
procedures incorporated in this analysis are: 

(1) The cloud location and volume are directly interpolated from the 
photographs of the stabilized cloud taken by three Askania tracking cameras. 

An ellipsoid shape for the cloud in each sector is assumed, 

(2) A Gaussian concentration distribution is used for each ellipsoid 
cloud. The concentration at the edge of the cloud is assumed to be 1% of 
that at the center. 
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(3) The total amount of pollutant is taken to be that resulting from 
17 sec of exhaust. The number 17 was chosen since early airborne samples 
show that the average HC1 monitored is about 4 ppm; this concentration, when 
combined with the observed volume of the cloud, 2.2 x 10 9 m 3 , gives a total 
mass of HC1 in the stabilized cloud of about 15 x 10 6 g , which equals that 
from about 17 sec of exhaust. In addition, from the camera observations, 
the cloud left in the PBL is that portion of exhaust which reaches the 
ground between firing and the attainment of 500 m altitude. Using the tra- 
jectory function for Titan III, 17 sec is required for the rocket to reach 
the altitude 500 m. 

(4) Since the camera observations show that the cloud contains two 
cloud portions, amounts corresponding to 14 sec and 3 sec of exhaust are 
assigned to the higher and lower portions respectively; this distribution 
ratio is based solely on the observed volumes of the two portions. 

Figure 26 shows the comparative results of the maximum ground dosages 
obtained using MDM and DISF. Compared to the results of calculations using 
the modeled cloud as input (Fig. 14) the results based on actual observed 
cloud characteristics show lower peaks and more rapid convergence of predic- 
tions from MDM and DISF at downwind distances. A comparison between this set 
of DISF calculations and airborne measurements is shown in Fig. 27 and a com- 
parison of the ground measurements and model predictions at monitoring site 
P-10 is shown in Table XI. Better agreement between airborne measurements 
and calculated results using the observed cloud rather than the modeled cloud 
(in Figs. 27 and 23, respectively) justifies our efforts. The maximum ground 
concentration path predicted from DISF using improved observed cloud charac- 
teristics as input is close to the ADPIC predictions in Fig. 15; we have pre- 
sented the predicted results before and after adjustment to the MDM predicted 
path as described in the previous section on ground-based measurement 
comparison. 
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F. DISCUSSION 


The comparisons made during this program demonstrate that the MDM general- 
ly overpredicts concentrations (or dosages) of exhaust constituents at the 
ground. This occurs principally because the model assumes that the standard 
deviation of wind azimuth angle a ^ is vertically uniform, whereas it varies 
with height. The concentration predictions are about an order of magnitude 
higher than the values indicated by the use of other models which give results 
in better agreement with the available data. For example, the DISF model gave 
good agreement with the data for the daylight August launch. The difference 
in concentrations predicted by MDM and DISF is about an order of magnitude, the 
same difference between TREATS and MDM for the December nighttime launch case, 
where TREATS gave values of the same order of magnitude as the measured data. 
Therefore, the ground level pollutant concentration predictions from MDM may be 
treated as conservatively high. 

As shown, the input of a is a dominant parameter for the MDM predictions. 
Direct measurements of vary rapidly over wide ranges, especially in night- 
time cases; it appears more useful to utilize the value of a ^ given by an in- 
terpolation method based on similarity theory which gives consistently higher 
predictions for a variety of meteorological conditions. 

There is no detailed information presently available on the cloud at sta- 
bilization; the description of the stabilized cloud given by MDM (more pre- 
cisely, by its preprocessor) seems to be reasonable when the center of the 
cloud mass is below the top of the PBL, although even under these conditions, 
MDM still tends to overpredict ground concentrations and dosages (see Section 
IV. E) . On the other hand, when the cloud center penetrates the top of the PBL, 
which usually happens with a shallow PBL (s 700 m) the total calculated source 
strength seems to be lower than that in the real cloud. This is a serious 
problem, because in this case the pollutant level is artificially decreased 
and the assumption of conservatively high estimates may be nullified. 

Recent efforts by SAI to modify the preprocessor program to more exactly 
account for the mixing rate of the exhaust plume and the rocket trajectory 
resulted in even greater overpredictions for ground concentrations. 

The following conclusions were reached concerning the other models used 
in this study: 
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(1) The moment scheme (TREATS) , coupled with the Blackadar formula for 
the vertical diffuslvity and the horizontal eddy dif fusivities de- 
duced from the Pasquill-Giff ord curves, did not properly model the 
cloud diffusion in the August case, but gave favorable comparisons 
with data for the nighttime December case. The probable explanation 
is that the eddy diff usivities used may not adequately describe the 
deep PBL (August case) which is driven mainly by sea breeze rather 
than the usual thermal and mechanical forces. The Blackadar formula 
and Pasquill-Gif ford curves having been based on the lower portion of 
the PBL, would be inadequate for the deep PBL existing during the 
August launch but valid for the shallow PBL of the December case. 

(2) The inherent ability of ADPIC to treat the advection effects due to 
wind variation is clearly exhibited in this study, particularly in 
the December nighttime case. However, this model showed very slow 
diffusion of the exhaust clouds; this is probably due to its implicit 
assumption of incompressibility of the pseudovelocity flow field 
which ignores cell expansion and consequently reduces the diffusion 
rate. However, some of the Lawrence Livermore Laboratory researchers 
believe that the grid was too large to resolve the source well — the 
numerical technique in these cases acted to further dampen the diffu- 
sion. 

(3) DISF gives the best agreement with the data for the August daylight 
launch case probably because of its use of off-diagonal dif fusivities 
and its physically sound Lagrangian approach to the diffusion param- 
eters. However, the predictions of this model do not correlate very 
well with data for the December nighttime launch case, especially in 
comparisons of the path of maximum ground level concentrations. The 
probable cause is the inability of DISF to physically describe wind 
direction changes. This is not an important factor in the August 
case with its low level of advection wind and high turbulence but it 
is important in the December case which was characterized by high 
wind and low turbulence. 
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( 4 ) 


The limited calculational results from METS do not provide sufficient 
information to judge its capability. However, its extremely low pre- 
dictions for the dosage field on the ground for the December case in- 
dicate the need for a more detailed examination of this model. 



V, DRY DEPOSITION AND AEROSOL COAGULATION STUDY 


During its initial history the Space Shuttle ground cloud contains 
both an aerosol, chiefly alumina, emitted by the solid state rocket motors, 
and a large amount of debris swept up by the rocket exhaust. In order to 
describe the evolution of particulate matter in the cloud during its early 
history (before cooling and condensation occurs), the following model of 
particle growth via coagulation including (i) agglomeration among the exhaust 
alumina particles, (ii) agglomeration between debris and alumina, and 
(iii) agglomeration among the debris particles has been formulated. The 
model allows the chemical composition of the evolving alumina/debris aerosol 
(debris is considered to be a homogeneous substance) and the rate of sedimenta- 
tion of this aerosol to be followed, 


A , FORMULATION OF THE PROBLEM 


Consider an aerosol containing (i) a chemical component made up of 
particles containing both a material A and a material B, referred to as com- 
ponent A/B, and (ii) a component of pure B particles, called B. 

Both components are divided into sets of particles having discrete 
masses, m^, where 


m^ = mi i = 1,2, ...q (46) 

where a is the mass fraction ratio. The number densities (particles ml -1 ) of 

A/ B 

the two components are then described by the sets of numbers {n^ (m^)} for 

component A/B and {n^ (m^)} for component B. Aerosol A/B is also described by a 
set of numbers (b^} where b^ is the average mass fraction of material B in ith- 
sized particles of A/B (0 < b^ < 1) . In order to obtain the set of particle 
number densities (n^(m^)} from a general particle number distribution, n(m) , 
one uses 


r (m i+i m i )1/2 ... f m ± al/2 

n- = I . n(m) dm = I n 

1 m .^ m .) 1 / 2 ^ a - 1 / 2 


i(m)dm 


(47) 


* When discussing properties of distributions, which are generally true for both 
A/B and B components, no superscripts are used. 
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If, as is often the case, n(m) is described by a power law 


n(m) = n 0 m' 


-x 


(particles ml” 1 g“; 


then 


where 


n^ = n(m£)m^ 


x 


sinh(x r log a) 


x-1 


The mass density of particles in the ith mass class is given by 


l/2 

£ i a ' C m i a 

n(m)mdm = n^ + / 

x a~ 7 •'m^a 


i« l/: 


n(m)(m-m-L) dm 


or, for the power law particle distribution. 


(48) 

(49) 

(50) 


(51) 


C m.a 

m . /*v 


1/2 




n(m)mdm 




1/2 


sinh[ (x’-l/2)log a] 
sinh[x'log a] 


(52) 


For rough estimates of the mass density of ith-sized particles, the term in 
brackets in Eq. (52) is nearly unity and the product may be used; for 

a % 10 and x in the range of 1.0 to 2.0 errors of up to 30% may occur if 
n-jm^ is assumed to be the particle mass in the ith mass class. Since the 
calculations which describe the coagulation process (see Eq. (74) below) 
rigorously conserve mass (if sedimentation losses are negligible) , any error 
in mass introduced by dividing the (continuum) initial mass spectrum into 
discrete mass classes is maintained throughout the calculation. If the 
shape of a power law distribution does not change significantly, the correct 
mass in each class can be recovered accurately from Eq. (52). Generally, 
however, one must reconstruct a continuum distribution n(m) from the n^ and 
then use Eq. (51) to obtain the actual mass in each class. 

In order to obtain expressions for the rate of change of the 
particle density, the following rates are defined: 
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L. . 
ij 


rate of collision between particles of component A/B, size i, 
with particles of A/B, size j. = twice the actual colli- 

sion rate.) 



I 


L' 


ij 


L" 


ij 


Collision 


= rate of collision of particles of component B, size i, with 
particles of component B, size j. (L' = twice the actual 

collision rate. ) 

= rate of collision of particles of component A/B, size i, with 
particles of component B, size j. 
kernels, K(m^, mj) = K , are defined by 


L. . 
13 

= K?. 

13 

A/B A/B 
n. n. 

1 3 

(53) 

L’ij 

= K? . 
13 

B B 
n. n. 

1 3 

(54) 

L". . 
13 

= K*7 . 

13 

A/B B 
n . n . 

1 3 

(55) 

form of 

the collision kernels has been given by Fuchs. 52 

The collision 


kernels considered here incorporate the effects of two processes: (i) the 

rate of Brownian coagulation and (ii) the rate of scavenging of light parti- 
cles by falling heavy particles, i.e. , 


K'. 

13 


= V 


13 


+ K, 


S.ij 


(56) 


K c , the contribution to the kernel from heavy particle sedimentation, is 
6 » ij 

given by 


K 


S,ij 


it C . . r v_ . - v c . 
13 max 1 S,i S,j 


(57) 


where r max is the radius of the larger of the two particle sizes 
where 



(r ± or r^ ) , 
(58) 


and C . . is a factor accounting for the fact that small particles slip around 
falling larger particles. 52 



(1 


+ r . /r )‘ 
mm max 


1 

1 + r . /r 

mm max 


(59) 


where r . is the smaller of r. or r.. The particle density (assumed the 
mm 13 

same for all particles) is denoted by p; the sedimentation velocity v s ^ of 
small particles is given by 
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*i' D l 


V S,i kgT (60} 

where is the particle diffusion coefficient (see below), k T , is Boltzmann's 
constant (1.38 x 10” 16 ergs K -1 ) and T is temperature. For larger falling 
particles, for which the particle Reynolds number approaches and exceeds 
unity, the boundary flow around the particle separates and drag increases 
more rapidly with particle size than is implied by Eq. (60). In this case 
the sedimentation velocity is obtained from an empirical relation given by 
Fuchs 53 : 



24 

Re^ 


+ 



(61) 


where the drag coefficient C 


D,i 


is given by 



2 


* Pg 




(62) 


and the particle Reynolds number is 

Re . = 

l 


2 r.p v. . 
i g 


(63) 


Using Eqs. (62) and (63) and Eqs. (64) and (65) which follow, for the 
density and viscosity of air (Pg and y, respectively) 

p g = 0.353 (P/T) (g ml" 1 ) (64) 

y = 1.055 x 10 -5 T 1 / 2 (g cm sec -1 ) (65) 

where P is in atm and T in K, Eq. (61) may be written in terms of P, T, v ., 

S , l 

r. , and m. as 

l i 

v_ . + 2.74 x 10 2 P 2 / 3 r. 2 / 3 T -1 v c . s / 3 - 4.93 x 10 s m,- r. -1 T~ l/2 = 0 
S,i l S , i 1 l 

( 66 ) 

In practice, the simpler Eq. (60) is used for particles with r^ <_ 25 ym and 

Eq. (66) is solved to obtain v . for larger particles. 

b , l 

The Brownian coagulation kernel is given by 2 
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r 


r. + r. 


4 (D. + D.) 

W m j ) ‘ 4lr(r l + r ]> < D i + V * + V A, . + (r. $ r ) 3 g. . 

' i J ij i m xj 


( 67 ) 


where, the diffusion coefficient. 


D = ^ 

x 6 tt yr 


- |^ 1 + 1.257 Kn ± + 0.40 Kn ± exp ^ J 


(68) 


and 


Kn . 


_ gas molecule mean free path _ 1 


particle radius 


r . 
x 


G„ is the mean thermal speed of ith and jth sized particles, 

G.. = {G 2 + G. 2 } 1 / 2 

xj i J 


(69) 


/ 8 kT X 1 / 2 
y i { it m i / 


(70) 


<5„ is a "mean free path" for particles: 


6 .. = ( 6. 2 + S . 2 ) 1 / 2 

iJ i J 


(71) 


(2) 1 / 2 tt G_. 
6 i = 48 D 

x x 


G. , / 8 D. . 3 r 8 D. . 2 : 3/ \ 

e{K + For) -K 2 + (fc7)] 


(72) 


For air at a temperature T(K) and pressure P(atm), the expression for the 
molecular mean free path, X, is: 


X = 2.2 x 10- 8 (|) cm (73) 

In order that mass be conserved and that particle number be 

5 £* 

accurately counted, the following method due to Kritz is used to divide 
the mass of a newly formed particle among the mass classes: If a particle 

of mass mj collides with one of mass m^ ( i j ) , the resultant particle is 

divided and a fraction F „ assigned to mass class i and a fraction (1 - F.jj) 
assigned to mass class i + 1. Since mass must be conserved 
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F..m. + (1 - = m. + m. 

13 1 13 l+l 1 3 


( 74 ) 


or 


m - — (m. + m. ) 

F . . - - i 1- - 1- 

Ij ” 1+1 - 


m. 


(a - l)m. 


(75) 


Using the expression for collision rates developed above, changes which occur 

a/b b 

in {n." /0 }, {n.°} and {b.l in a time interval At may be calculated. Changes 
due to coagulation are computed as follows (losses due to sedimentation from 
the cloud will be taken into account later) : 


A /B - f A 


1 

(L.. + L"..) + 2(S..L.. + L"..) F.. 

13 13 3=1 13 13 1 3 13 


i-1 

+ 2 [(S. n ,L. . . + L". ,)(1 - F. . .) + L" . . F. . ] 

. =1 i-1, 3 i-l,J i"l,3 1-1,3 31 13 


i- 2 

A A, 1-1 (1 -wj 


At 


( 76 ) 


A n 


B 


{ 


i-1 i 

- 2 (L* . . + L" . . ) + 2 S. .L\ .F. . 

j-l 1J 31 j-l 13 13 13 


i-1 


+2 S. . .L' . (1 - F 

j =1 i-l,3 i-l,3 


i-l,j } j 


At 


( 77 ) 


where 


ij 


c 


1/2 i=j 
i#3 


( 78 ) 


The factor S.. is required for correct counting. L. . and L' . . are 

13 11 11 

twice the actual collision rate for ith-sized particles with one another. 

Each such collision removes two of these particles from the ith class and 

places a fraction F. . of a particle in the ith class and (1 - F..) in the 

IX xx 

(i + l)th class. 
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The mass of material B added to the ith class of A/B in time At is 


given by: 

am. a/b (b) 


{ 


q i 

- 2 (L.. + L"..) m.b. + 2 [S..L.. (m.b. + m.b.) 

iJ iJ 1 i y y 1 i 3 3 


m. 


+ L". . (m.b. + m.)] ^ F. . 

13 i i 3 m i + m j *3 


i-1 

+ 2 [s . L . (m. 1 b. + m.b.) 

j =1 i-l»3 1-1.3 i-l i-l 3 3 


+ L 


m # 

". .. . (m. ,b. + m.) ] 

1-1,3 i-l i-l 3 m ± _ 1 + 


(1 - 


i-l m. 

+ 2 LV . (m. + m.b.) F. . 

j=1 3 ,i i 3 3 + m 2.3 


i-2 \ 

+ 2 L n . .. (m. n + m.b.) m i n „ x l 
j-1 J,i-1 1-1 3 / — + — (1 - F HiJ ) > 


X At 


(79) 


To calculate the change in b^ in At, 


Ab. = 

l 


M. A/,B (B) + AM. A/B (B) M. A//B (B) 

1 1 _ _1 

, A/B . A/B. 
m . (n . ' + An . ) 

x i x 


m. n . 
x x 


A/B 


(80) 


or 


Ab 


{ 


1 + AM. A/B (B)/m.n. A/B b. 
x xx x 

, , A A/B, A/B 
1 + An. /n. 
x x 


h 


(81) 


a/b b a/B 

To obtain the total changes in the {n. } and {n. }, {An. } and 

{An_^ }, we must add to the changes caused by coagulation, those caused by 
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sedimentation. These latter losses, assuming the cloud to be well-stirred 
(homogeneous) , are given by 


A/B 

A"n A/ ^ B = — — S - ^ — — At 

i h 

B 

„ v„ . n. 

A" n . B = - At 

x h 

where h is the cloud height. The total changes are then given by 


(82) 


(83) 


a„. a/b 

X 

- A’n. A/B + A"„ A/B 

X X 

( 84 ) 

CQ 

•H 

£ 

At 8 , A || B 

= An. + An. 

X X 

( 85 ) 


It remains to specify At. In order for the calculation to accurately 

follow the evolution of the aerosol, An./n. must be much less than 1 for all 

x x 

i. It can be seen that (Eqs. (76), (77)) this will be the case if 

iK?. n. - K?. n. F. . i At < < 1/q, j ^ i 

l 13 1 i] 3 i]) 


( 86 ) 


and 


K. . n. At « 1/q 
ij 1 


3 > i 


for all i (and both B and A/B aerosol components). Thus one sets 

At = min(Ati, At 2 ) 

where 


(87) 


At x = 


Tq - 1 ( ^ Tq 

i-n.\ /d - F..) K? n\ 

l m i J Jmax 13 max 


for 


a ^ 2 j 


3 — 1 


(§ 8 ) 


where the subscript (max) indicates the maximum for all i (and all j < i ) . 
Also , 


. . _ Tq 

Atz {K9. n .} 

xj 3 max 


3 > i 


(89) 
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where (max) indicates the maximum for all i (and j > i) and the time step 
size parameter, T, is chosen in such a manner (i.e., large enough) that times 
required for computation are practical, but small enough that the desired 
accuracy is obtained. 


E, CODE VERIFICATION 

The equations developed in the previous section were coded and a 
series of calculations was performed on model systems in order to test how 
accurately the evolution of a cloud of coagulating particles is followed and 
what effects the time step and mass class size parameters, T and a, have on 
this accuracy. 

Figures 28 and 29 show the results of calculations, performed on a 
model, single component system for which the integro-dif f erential equation 
governing the rate of coagulation, 52 


3 

3t 


n(m, 


t) 



K n(m - m' , 


t) n(m' , 


t) dm' 


CO 

K n(m, t) n(m', t) dm’ (90) 

has an analytical solution. This model system is one in which K is a constant 
and the initial (t = 0) distribution is 



n(m, 0) 


rip m 
m 0 2 


exp (—2 m mo -1 ) 


(91) 


where n 0 is the initial (total) particle number density and m 0 is the initial 
average particle mass. At later times, the particle distribution follows the 
equation 


n(m, t) = 8 


n 0 exp (—2 m m 0 1 ) 

mo(noKt) 1 / 2 (n 0 Kt + 2 ) 3/2 


sinh 2 m m 0 


" ( < 92 > 


For the case considered in Figs. 28 and 29, K = 6 x IO -10 ml particles -1 sec -1 , 
n 0 = 1 x 10 6 particles ml -1 and m 0 = 5.24 x IO -10 g (corresponding to a 10 ym 
particle with unity density); times of t = 0, 1.25 x 10 s sec, and 1.7 x 10 6 
sec are considered which correspond to a dimensionless time (K n 0 t) of 0, 75 
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and 1020. The analytical solution of Eq. (90) in the form of 


N(d,t) = particles ml - with diam > d 


= f\ 


On) d(d) (93) 


and 


N’(d,t) = particles ml -1 with diam 


r d 

iam < d -I n 


(m) 


dm 

d(d) 


d(d) (94) 


are shown as curves in Figs. 28 and 29, respectively, 
machine calculations giving values for 

V (d i < W l/a) - \* 3 

3=i+l J 


These are compared to 


(95) 


and 


, s \ ^ 

J =1 


(96) 


which were performed using the values for the mass class ratio, a, and the 
time step parameter, T, indicated in the figure captions and also in Table XII 
which summarizes the results shown in Figs. 28 and 29. 

Figures 28 and 29 and Table XII show that the numerical solutions 
yield, to a good approximation, the actual solution to the coagulation rate 
equation. Numerical solutions were obtained using the following pairs of 
values for the mass class ratio, a, and the time step parameter, t: a = 2, 

T = 0.01; a = 2, T = 0.1 and a = 10, T = 0.1. The total number of particles 
and the mean mass or, equivalently the mass-average diameter d^ = 3 

are accurately calculated.* However, it is apparent from Figs. 28 and 29 that 


* The mean mass at time t is obtained from the output by the formula 


m(t) 


/ r . r 

(2 n. (t) m / 2 
' 3=1 i=1 


ru (t) 



M 


total 


r 


/.^i 

i=l 



i.e., it is assumed that the second term on the right hand side of Eq. (51 ) 
remains a constant percentage of the first term. Here M tota 2 is the actual 

initial mass of the aerosol as given by the left hand side of Eq. (51). 
Since, if no sedimentation losses occur, mass is conserved in the calcula- 
tions, the correct calculation of m(t) depends only on the correct calcula- 
tion of the total particle number. 
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the numerical solution tends to overestimate the number of light and very heavy 
particles. The tendency becomes very noticeable for the larger mass class 
ratio of a = 10. Since mass is rigorously conserved in the numerical calcula- 
tion, the overestimation of numbers of heavy particles necessitates that the 
mass of the remaining lighter particles be underestimated. This can be seen 
clearly in Fig. 29 and in Table XII where the values of d N qq, d N ^Q and d^iOj 
the particle diameters for which N (Eq. 93) is 90%, 50%, and 10% of the total 
number of particles, are all found to be smaller than their true values. 

For the user of this code interested in keeping computational time 
to a minimum it is important to note that the use of x = 0.1 instead of 
T = 0.01 with a = 2 leads to no noticeable decrease in accuracy. On the 
other hand, increasing a from 2 to 10 leads to a marked decrease in accuracy; 
with a = 10 the numerical solution (while still having accurate values for 
total particle number and dgj-) has become badly "smeared" with, e.g., d N 
the median particle diameter, being only 70% of the correct value at t = 

1.7 x 10 6 sec (see Table XII). Since calculations such as these can often 
require several minutes on a modern computer (a CDC 6600 in this case) it is 
clear that setting T as large as possible (x 0.1) is desirable since 
little or no accuracy is lost. On the other hand setting a much greater 
than 2 results in a very substantial loss in accuracy and one must then 
weigh this loss against the cost of machine time. 

Before suggesting what would appear to be good practical values 
for X and or it is worthwhile to examine the results of calculations on a less 
artificial system. In Fig. 30 the results of calculations on a system in 
which the coagulation of an aerosol initially containing 1 x 10 s particles 
ml -1 is followed. The initial size spectrum of the aerosol is of a form 
often obtained in practice, 

n(m) = n£ m -2 5.24 x 10 -19 g < m < 5.24 x 10~ 13 g 

or 

n(d) = n 0 d -z * 0.01 ym < d < 1.0 ym 

c 

In this case the collision kernel K. . is of the form given by Eqs. (56), (57) 
and (67) and, for selected values of particle masses, is tabulated in 
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Table XIII. No analytical solution of the coagulation rate equation (Eq. (90)) 
can be obtained for this case and we will assume that the numerical solution 
using a = 2.07 and T = 0.01 is accurate. As with the case just discussed, it 
can be seen that no great error is incurred by increasing T to 0.1 and indeed, 
in the calculation at t = 250 sec, no major errors (assuming here that the 
T =0.1 calculation is correct) appear when T is increased to 1.0 except for 
the two lightest mass classes. (At heavy masses, d > 1pm, there are no 
significant differences in any of the cases.) For a = 10 at 130 sec, errors 
are very apparent at light masses, but no significant error is noted in the 
total particle number of the median particle size (d^ ,-q) . At t = 1.3 x 10** 
sec only results for T = 1.0 are shown; excessive computational times would 
have been required to obtain results using T = 0.1. This curve, if the 
trends noted at t = 250 sec are maintained, underestimates the number of 
very small particles but should yield approximately correct values for total 
particle number and mass average particle size. 

From these examples the following "rules of thumb" regarding the 
selection of a and T are apparent. First, since values of T up to 1.0 have 
little effect on the sample calculations, it is probably true that use of 
values of 0.1 to 0.5 will not result in meaningful errors. Only where the 
size spectrum is very rapidly changing, as in the smallest size classes of 
Fig. 30 where small particles are rapidly being consumed, would any signifi- 
cant errors result. (The calculated rate of decrease of these very small 
particles has little or no physical significance anyway, since in any real 
system no sharp cut off at 0.01 ]im would exist.) Secondly, the size of a, 
the mass class size ratio, should be kept as small as possible consistent 
with economical computation. Use of a = 2 yields very accurate results. 

For a = 10, total number and mean particle mass are computed accurately but 
the shape of the calculated particle distribution may tend to become 
incorrect if a rapidly varying size spectrum such as the exponentially vary- 
ing one of Figs. 28 and 29 is treated. 

As a rule, then, the use of a in the range of 2 to at most 5 is 
indicated. It should be kept in mind that the computational time 
will be inversely proportional to a 2 and thus, if extended calculations are 
required, every effort should be made to make a as large as is consistent 
with the demands to be made on the results. A particular example of some 


70 



importance is that in which it is desired to accurately predict sedimentation 
losses. In this case the possible overprediction of the production of very 
large particles, such as was seen in the comparison with the analytical 
solution (Figs. 28 and 29), may lead to an overestimation of the amount of 
material lost. Unfortunately, time did not permit us to systematically quan- 
tify this observation; a few calculations where this comparison could be made 
indicated that sedimentation losses could be overestimated by as much as a 
factor of 2 when a = 10 was used instead of a = 2 with power law distributions 
of the type used in Fig, 30. 

Both of the verification tests shown in Figs. 28, 29 and 30 deal only 
with coagulation of a single component aerosol. In order to check that the 
mixing of two components is correctly treated the following simple test was 
conducted. A calculation was performed in which the initial aerosol was one 
containing two components, A and B, with identical size spectrum. As this 
aerosol coagulates, those particles which are mixtures of A and B should 
display equal amounts of material A and material B, i.e, the following rela- 
tionship should hold 


. n ,, A/B B w A/B 

b± = 0.5 (n ± — n i )/n i 


Calculations performed over several thousand time steps showed this relation- 
ship to hold rigorously. 


c • DRY DEPOSITION FROM THE GROUND CLOU D 

The coagulation code has been used to compute the loss of alumina 
from the ground cloud of the May 20 Titan launch during the cloud rise period. 
It is assumed that no particle growth by means other than coagulation occurs, 
i.e., the warm cloud remains above the dew point. The cloud is considered 
to be well-mixed and thus to have a uniform density of particulate matter 
throughout. The expansion of the cloud via entrainment and the pressure 
decrease during cloud rise, as discussed in Section II, are incorporated into 
the calculations. 

The calculations, performed over a time interval in which the cloud 
grows from an initial volume of 6.5 x IQ 7 m 3 to 5 x 10 9 m 3 when stabilization 
occurs at 525 sec at an altitude of about 1800 m, are obtained by using Briggs' 
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formula. The mass of alumina deposited in the cloud by the Titan is 2.2 x 
10 7 g and thus the initial density of particulate matter is 0.340 mg m -3 . The 
density of the alumina (and the swept-up debris) is taken to be 3.0 g ml -1 . 

Uncertainty as to the actual size spectra of the alumina and debris 
particles and the mass of debris swept up requires that a range of values for 
these parameters be treated. Calculations were thus performed (i) using three 
alumina distributions from different sources, (ii) assuming that, for particles 
less than 1 ym in size, the mass of the debris is 10 times the mass of the 
alumina, and (iii) assuming the total mass of debris is either 100 or 1000 
times that of the alumina. The three alumina particle size distributions are 

S 5 

those of Varsi, obtained from Titan flight measurements, and those of 
Dawborn 56 and Kreautle 57 which were obtained from measurements on particles 
produced in small motor firings in tanks. These distributions are given in 
Table XIV. The measurements of Kreautle include rapidly settling particles 
collected from the floor of the tank and are thus noticeably different from 
the other two measurements of airborne particles in having considerably more 
large ( ^ 1 ym) particles. Thus the Kreautle data would represent a 
worst case situation with regard to sedimentation early in the cloud history 
if no debris were present. (As will be seen below the presence of debris 
modifies this conclusion somewhat.) The assumption that, below 1 ym, there 
will be about 10 times as much debris as alumina is based on airborne measure- 
ments of the Titan ground cloud. The assumption that either 100 or 1000 times 
as much debris as alumina is to be found in the cloud initially is simply a 
guess.. The results of calculations presented here are for a factor of 1000 
for calculations with the Varsi and Dawborn alumina distributions and both 
100 and 1000 with the Kreautle alumina distribution. (It will be seen that 
the amount of alumina carried down is not overly sensitive to this ratio and 
thus the accuracy of these calculations is probably not limited by this 
unknown factor so much as by such simplifying assumptions and errors as that 
of assuming a uniform distribution of particulate material in the cloud.) 

Table XIV also summarizes the size distributions used for the debris. The 

median particle diameter d„ T and the mass-average diameter d_ of all 
r N,50 m 

distributions are also given. The calculations were performed using a time 
step parameter, T = 0.01, and a mass class size ratio, a = 10. As discussed 
above, this results in some overestimation of the rate of formation of large 
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particles and thus tends to be a worst case calculation with regard to the 
amount of sedimentation. Figures 31 to 34 show the particle number 
distribution, N r (d) (particle with diam < d) and mass distribution, M(d) (mass 
of particles with diam < d) for total mass and alumina for the four test calcu- 
lations. Note that for the Kreautle distribution the initial mass distribu- 
tion for alumina is very nearly a step function, dropping rapidly from its 
maximum at 340 mg m -3 to zero between 2 and 7 ym. The Varsi and Dawborn ini- 
tial distributions have more curvature with most of the mass being contained 
in particles less than 1 ym in size. 

At stabilization the cloud has expanded by a factor of 72. Total 
particle number densities have dropped factors of 250 to 600. As can also be 
seen in Figs. 31 to 34 the larger debris particles which are still airborne 
at stabilization have captured about 2% to 5% of the alumina. This effect is 
denoted by the rather abrupt change in slope of the M curves occurring at 
d ^5 ym in the figures. (In the Varsi data case the change is not so 
noticeable as in the other three cases.) This process of scavenging of small 
alumina particles by large, rapidly falling, debris particles is responsible 
for most of the removal of alumina from the cloud. For example, using the 
Kreautle alumina distribution with the debris mass being 1000 times the alu- 
mina mass at 525 sec (stabilization) the rate at which alumina is carried 
down from the cloud by large debris particles with diameter > 10 ym is about 
10^ times the rate of loss of alumina by the sedimentation of pure alumina 
or small alumina/debris particles. It is important to note from these figures 
that the amount of alumina coating the larger debris particles is not very 
dependent (i) on the particular alumina distribution used or (ii) on the 
amount of debris present. This is further exemplified in Figs. 35 and 36. 

In these figures the amounts of sedimentation, both total and of alumina, 
are plotted. The Kreautle distribution with 1000 times as much debris as 
alumina and the Dawborn distribution with the same mass of debris lead to 
the deposition of 3 x 10 6 g and 2 x 10 6 g, respectively, of alumina in the 
first 525 sec (375 sec for the Dawborn case). The use of the Varsi distribu- 
tion with 1000 times as much debris and Kreautle distribution with 100 times 
as much debris as alumina, deposit about 1 x 10 6 g. The greater deposition 
from the Kreautle distribution with the larger amount of debris is expected; 
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thus that the Dawborn distribution would yield about the same amount of 
deposited alumina is not so easy to explain. Apparently the amount of 
material which settles is a relatively strong function of the average size 
of the debris particle (more so, for example, than of debris particle number) 
which is larger for the Dawborn and Kreautle (1000) distributions. In any 
case, the differences are not great between any of the cases. However, if 
debris had not been present, large differences would have been noted with 
the Kreautle distribution yielding about an order of magnitude more sedimenta- 
tion than the other two. 

If one assumes that the stabilized cloud were stationary, spherical 
and uniform with a volume of 5 x 10 9 m 3 and the deposited alumina were spread 
evenly over the ground under the cloud, the worst case loading of alumina, 
averaged over the ground beneath the cloud would be 0.8 g m -2 at the end 
of 525 sec. The drift of the cloud will, of course, tend to decrease actual 
loadings , but the nonhomogeneity of the cloud and the fact that early in 
the cloud's history it is smaller would tend to raise the loading level near 
the launch pad. Although the calculations do not take into account the 
movement of the cloud, which tends to reduce loading at any particular point 
on the ground, the rate of loss of material due to sedimentation throughout 
the cloud rise period is calculated. It is possible to integrate this rate 
and thus calculate the loading at the point on the ground below the center 
of the (stationary) cloud. These maximum time-integrated loadings, at the 
end of the cloud rise period, are given in Table XV. The trends are similar 
to those for total loadings of Figs. 35 and 36. The Kreautle distributions 
with heavy debris loading give the heaviest alumina loading at the ground, 

5 times that of the Kreautle distribution with less debris. 
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VI . CONCLU SI ON S AND RECOMMENDATIONS 
A ■ CONC LUSI ONS 

This study, which was undertaken to investigate the complex interactions 
between the turbulent flow field in the lower portion of the troposphere (the 
PBL) and the rocket exhaust ground cloud resulting from Space Shuttle/Titan III 
rocket launches, comprises two separate tasks: (1) an attempt to identify and 

minimize the uncertainties and potential inaccuracies of the NASA Multilayer 
Diffusion Models using data from selected Titan III launches at KSC, and (2) a 
systematic analysis of the physical/chemical processes taking place during the 
cloud rise, and formulation of a realistic time-dependent model. The former 
study is based on detailed parametric calculations using the MDM code and a 
comparative study of several more exactly formulated diffusion models, the 
MDM, and NASA measurements. 

The results of the comparative studies of diffusion models and the para- 
metric calculations are as follows: 

(1) If the input standard deviation of the azimuth angle, o^, is chosen 
appropriately, the MDM consistently overpredicts the ground level 
concentrations and dosages for the cases examined in this study (cf. 
data and other models). However, it should be noted that the MDM 
predicted results at ground level are strong functions of a ; the 

A 

peak maximum ground level concentration (or dosage) predicted by MDM, 
using an interpolated o based on similarity theory and assuming that 

A 

is vertically uniform in the entire PBL, is an order of magnitude 
larger than the value indicated by the available data. Therefore the 
peak in maximum ground level concentrations predicted by MDM seems to 
be a conservative input for analyzing ground level environmental haz- 
ards (such as the environmental constraint for Shuttle launches) . 

(2) MDM overpredictions at ground level are mainly due to rapid transport 
of more highly concentrated pollutants from the upper portion of the 
PBL. This results from the assumption of strong turbulence (large 
values of a ) in the upper portion of the PBL. Conversely, the pos- 

A 

sibility must therefore be considered that pollutant concentrations 
above the ground level may be underestimated. Environmental hazard 
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analyses, e.g. weather modification studies, which require, as input, 
pollutant concentrations throughout the entire PBL, must employ MDM 
with caution. 

(3) In a shallow PBL condition, the present study shows that the uncer- 
tainty in the entrainment constant used in the cloud rise formula 
can cause a factor of three (2 3) error in the downwind ground level 
concentration predictions (Section IIL A. 2) . It is also found under 
such PBL conditions that the MDM-simulated initial cloud contains 
lower pollutant concentrations than does the real cloud. Although 
the latter deficiency may occasionally nullify the "conservative" 
prediction ((1), above) the problem can be resolved by simply assum- 
ing that the center of the exhaust cloud mass (stabilization height) 
is located below the PBL height. 

(4) Compared to the more exactly formulated models, TREATS, ADPIC, and 
DISF, the MDM model has less potential for including nonhomogeneous 
and nonstationary features such as variations of wind direction and 
wind speed, non-uniform surfaces, land-sea interactions, etc., or for 
incorporating microphysical and chemical processes (such as aerosol 
mechanics and heterogeneous chemistry as well as kinetic chemistry). 
Nevertheless, it should be kept in mind that the MDM model, because 
of its ability to consistently overpredict ground level concentra- 
tions, can certainly serve as an acceptable engineering tool for use 
in environmental hazard analyses at ground level. 

(5) The TREATS, ADPIC and DISF models although they serve only as refer- 
ence tools in this study, do contribute to an understanding of the 
relative merits of the different diffusion modeling techniques. For 
example, the moment scheme of TREATS demonstrates the mathematical 
relationship between the vertical nonhomogeneous nature of turbu- 
lence and the moments of the pollutant concentration distribution; 
the trajectory technique of ADPIC provides a description of advec- 
tion effects due to wind variations in time and space. Although the 
model in its present form predicts lower pollutant concentrations at 
ground level than indicated by available data, it can give a good 
description of the exhaust cloud transit path. The ADPIC-predicted 



ground level pollutant concentration can also be treated as a lower 
bound. The DISF model demonstrated the necessity of including off- 
diagonal diffusivities and indicated the advantages of the physically 
sound Lagrangian approach for the derivation of diffusion parameters. 
The overall merit of the DISF model has been shown in its favorable 
comparisons with airborne measurements in a strong turbulence condi- 
tion. The various strengths of these three models may serve as valu- 
able guidelines for the future development of a new and advanced dif- 
fusion model. 

(6) The current lack of micrometeorological information concerning the PBL 
at KSC causes many difficulties in making diffusion calculations and 
results in wide uncertainties in the calculated results. Based on the 
MDM parametric study, the MDM-predicted results are found to be more 
sensitive to those parameters relating to the micrometeorological in- 
formation, such as o and PBL height than to uncertainties in the de- 
scription of initial cloud and ground surface effects. For example, 
the peak of the maximum ground level concentration from MDM can vary 
by an order of magnitude depending on the value of a (Fig. 4) 
chosen. A factor of two error is encountered due to the uncertain- 
ties of the PBL height. 

The major conclusions of the calculations on the physical/chemical pro- 
cesses during cloud rise are as follows: 

(1) The value used in the MDM for the heat release from the Shuttle solid 
fuel (or Titan III) is reasonable. 

(2) Deluge water injected into the exhaust plume from the AWSS for the 
Shuttle, or from the trench for Titan, has little effect on subse- 
quent concentration predictions if wet chemistry is not an important 
factor. 

— 2 

(3) It was found that an average loading of about 1 g m of alumina 
(about 15% of the alumina initially in the cloud) , with respect to 
the particle size distributions tested, will be deposited on the 
ground in a worst case example, as in the Titan III launch of 20 May 
1975. 
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(4) Beneath the center of a rising cloud with no horizontal movement, load 
ings as large as 4.5 g m -2 are calculated. Horizontal movement of the 
cloud will lower the loading and thus there is good reason to believe 
that the calculated results include the worst case situation. 

(5) In dry deposition alumina settles predominantly via coagulation with 
rapidly falling debris particles. While more information on the 
alumina and debris particle size distributions and the mass of debris 
would be useful, the amount of alumina deposition is not critically 
sensitive to these parameters. 

B . RECOMMENDATIONS 

The following recommendations are based on our study of the interactions 
between the rocket ground cloud and the atmosphere: 

(1) This comparative study of diffusion models focused on a sea breeze 
daylight case and a nighttime case. An extension of this type of 
study to additional meteorological conditions would be useful. 

(2) At this stage, the models adopted in this program, with the exception 
of the MDM, have not been fully documented and some of them (e.g. 
ADPIC) can only be used in the developer's computer system. There- 
fore efforts should be made to fully document those models and/or to 
convert them into a more widely used computer system. 

(3) A simulation model to provide micrometeorological information on the 
PBL at KSC, with special emphasis on determining the diffusion para- 
meters and the thickness of the PBL, should be developed. This model 
must take into account the land-sea interaction feature at KSC, but 
should retain its less complex features for routine usage. 

(4) This study showed that the calculation and comparisons which use, as 
input, a description of the stabilized cloud based on observed cloud 
location and volume, gave better agreement with the data than those 
which use the MDM modeled cloud. Therefore, an engineering formula 
or modifications of the cloud model in MDM on the basis of empirical 
correlations between the meteorological characteristics and available 
observed data of Titan III clouds should be implemented. 
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(5) Heat flux, v'q^ near ground height (say 4 m) should be measured at 
existing tower locations at KSC, because it is a fundamental quan- 
tity in determining the physical scales in the PBL, such as charac- 
teristic wind speed, characteristic temperature, Monin-Obukhov 
length, etc. 

(6) At present, the measurement data on the size distribution and com- 
position of aerosols associated with the ground cloud for validation 
and refinement of the models are scarce. Efforts to obtain such data 
should be encouraged. 

(7) In this study, only limited calculations on cloud dry deposition were 
made. Extended calculations for the various meteorological condi- 
tions and for different rocket exhausts such as from the Shuttle 
should be carried out. Some modifications and refinements of this 
model, accounting for the different densities of two kinds of parti- 
cles and turbulence- induced coagulation, and incorporating wet chem- 
istry (condensation) should be made. 
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APPENDIX A 


CLOUD RISE AEROSOL MODEL (CRAM) 

COMPUTER PROGRAM 

I . PROGRAM DESCRIPTION 

This program is used to calculate coagulation between two types of aerosols 
in a spherical buoyant atmospheric cloud, or in a group of buoyant clouds in 
which the lower clouds are spherical and the upper clouds are cylindrical. The 
changes in the cloud, i.e., the increase in size, and the decrease in tempera- 
ture and pressure, are taken into account. The aerosols are divided into a set 
number of discrete mass sizes. 

For calculation purposes, one aerosol is kept pure, while the other con- 
tains a fraction of materials from each type of aerosol. A separate fraction 
is kept for each mass class. Both sedimentation and Brownian coagulations are 
calculated, as well as sedimentation loss through fallout. If the multi-cloud 
option is used, sedimentation loss in one cloud is added to the cloud below. 

For generalization of application, two other options are included: 

(1) PTVC option: Pressure, temperature and volume are kept constant. 

(2) KKC option: The collision kernel is kept constant. 

Output provided historically: 

(1) Mass and size distribution for each aerosol. 

(2) Sedimentation loss: Sedimentation velocity, total mass loss, number 

lost in each mass class, and mass loss rate are included. 

(3) Physical properties of each cloud, including the entrained air mass 
from which the water content in the cloud can be calculated . 

II., PROGRAM INPUT 

A. FORMAT DESCRIPTION 

Input must be in the form: 

<title card> 

<namelist JIM> 

<namelist JIM1> 
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Where the title card is one line of alphanumeric data. Columns 2-80 are 
printed at the top of the output. 

Namelist data has the form: 

$<namelist name> <variable lists> $END 
Column 1 must always be left blank. The. variable list has the form: 
<variable name> = <constant>, or 

<array name> = <string of constants (separated by commas) > 

Every variable list (including the last one) must be followed by a comma. 
If two or more constants in a string have the same value, they can be repre- 
sented as n* (value), where n is the (integer) number of repetitions. Blanks 
must not appear inside variable or namelist names, but they can be used any- 
where else. Constant strings can be continued from one line to the next. 
Variables can be defined in any order inside the namelist, but the two name- 
lists must be in the proper order. Namelist JIM1 is not used if the PTVC 
option (constant pressure, temperature, and volume) is selected. 

If the variable has a default value and is not defined in the namelist, 
the default value is used by the program. 

If the multicase option is used, then complete sets of data, including 
the title card, must be given for each case. However, after the first set of 
data, input values do not have to be defined if they are the same as the 
values defined in the first set, except for the arrays ANP, BNP, and BNEW 
(which must be defined every time) and variables with a default value (which 
must be defined any time a value other than the default value is desired). 
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Value 


Array 


Default 


Name 

Type 

Units 

Restriction 

Size 

Value 




NAMELIST JIM 



ALFA 

real 



>2 




ANP 

real 

no-ml 1 

— 

30 x 20 

— 

BNEW 

real 

— 

— 

30 x 20 

— 

BNP 

real 

no-ml” 

— 

30 x 20 

— 

CNH 

real 

cm 

— 

— 

— 

CNOMX 

real 

no-ml 

— 

— 

— 

CNPRES 

real 

atm 

— 

— 

— 

CNTEM 

real 

°K 

— 

— 

— 

CNVOL 

real 

m 3 



„ . _ r . 



FACTRA 

real 

— 

— 

20 

1.0 

FACTRB 

real 

— 

— 

20 

1.0 

IOINT 

Integer 

— 

>1 

— 

— 

IR 

integer 

— 

<30 

— 

— 

KKC 

logical 

— 

— 

— 

.FALSE. 

NCLTI 

integer 

— 

— 

— 

— 

NLAY 

integer 

— 

3. £ n 5 20 

— 

1 

NT 

integer 

— 

— 

— 

— 

NT1 

integer 

— 

— 

— 

1000 

NT 2 

integer 

— 

— 

— 

1000 

NZS 

integer 

— 

— 

— 

1 

PTVC 

logical 

— 

— 

— 

.FALSE. 

RHO 

real 

g-cm” 

— 

— 

— 

RMO 

real 

g 

— 

— 

— 

RN2M 

real 

— 

— 

— 

1.0 

TEND 

real 

sec 

— 

— 

— 

TI 

real 

— 

— 

— 

— 

XKCON 

real 

ml-sec - 1 

— 

— 

— 


NAMELIST JIM1 


BZN 

real 

m 

CORTM 

logical 

— 

DPHZ 

real 

k-m” 1 

DT 

real 

sec 

DWATER 

real 

cal 

GAMTWO 

real 

— 

INTPR 

logical 

— 

IPINT 

integer 

— 

NLAY 

integer 

— 

P 

real 

atm 

QC 

real 

cal sec 

T 

real 

°K 



.FALSE. 


100 

100 


3.2557E9 

0.49 

.TRUE. 

5 

1 

1.1883E10 
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Value 


Array 


Name 

Type 

Units 

Restriction 

Size 

TF 

real 

sec 




TFA 

real 

— 

— 

— 

TFB 

real 

— 

— 

— 

U 

real 

m sec -1 

— 

100 

ZP , ZU , ZT 

real 

m 

— 

100 


Default 

Value 
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C. DESCRIPTION OF VARIABLES 


Name 


ALFA 

ANP , BNP , BNEW 

CNH 

CNOMX 

CNPRES 

CNTEM 

CNVOL 

FACTRA , FACTRB 

IOINT 

IR 

KKC 

NCLTI 

NLAY 

NT 
NT1 
NT 2 

NZS 

PTVC 

RHO 


Symbol Description 

NAMELIST JIM (MAIN PROGRAM) 


a 

B 
n . 

1 

A/B 
n . 


b i 


P 

T 


q 


p 


Mass ratio of consecutive mass sizes. 

Arrays giving A/B distributions, B distributions, 
and fraction B in A/B, respectively. Each array 
has 600 elements; 1-30 for the bottom cloud, 

31-60 for the second cloud, up to 20 clouds; 
elements 1, 31, 61, etc. give number for smallest 
mass particle in each cloud. 

Used only with PTVC option (normally calculated 
in CLDRIS routine). Vertical thickness of cloud. 
(Used only for calculating sedimentation loss). 
Sedimentation loss can be rendered insignifi- 
cant by making CNH very large. 

Used only with PTVC option. Minimum value for 
mass distributions. (Values less than CNOMX 
are set equal to CNOMX at the end of each time 
interval . ) 

Used only with PTVC option. Cloud atmospheric 
pressure. 

Used only with PTVC option. Cloud temperature. 
Used only with PTVC option. Cloud volume. 

At the start of the program, the A/B and B mass 
distributions of cloud number are multiplied by 
FACTRA (n) and FACTRB (n) , respectively, for all 
n up to NLAY. (This provides a convenient scal- 
ing factor.) 

Number of time intervals between printouts. 

Number of mass sizes. 

Logical variable. When TRUE, the collision 
kernel is set equal to a constant. 

Number of time intervals between recalculations 
of delta time ("DT"). 

Number of clouds (can be defined in CLDRIS 
routine, but it must be defined here if it is 
not equal to 1 and the PTVC option is used) . 
Number of time intervals to be calculated, if 
cloud rise time (TEND) is not exceeded first. 
Increment value for NT2 when time interval NT2 
is reached. 

Number of time intervals at which to first 
revise TI; then counter for subsequent revisions 
of TI. 

Multicase option. If NZS is not equal to 1, 

CRAM expects another complete case of input data. 
This is a logical variable. When it is TRUE, 
the constant pressure, volume, and temperature 
option is used. 

Density of particle material. 
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Name 


Symbol 


Description 


RMO 

RN2M 

TEND 

TI 

XKCON 


mi Mass of the smallest particle (smallest mass 

size) . 

Revised value for TI (TI=*RN2M) . 

Time at which to stop calculations (cloud rise 
time)* if NT not exceeded first, 
x Initial value for time step size index. 

K. . Value of collision kernel when KKC option is 

used. 


BZN 


CORTM 


DPHZ 

DT 


DWATER 

GAMTWO 

INTPR 

IPINT 

NLAY 

P,U,T 

QC 

TF 

TFA, TFB 
ZP, ZU, ZT 


NAMELIST JIM1 (CLDRIS ROUTINE) 


A<j>/Az 


Y 


P , U , 
amb , amb , 


lamb 


Q 

t. 


Array giving initial altitudes of bottom of 
each cloud, and top of top cloud, (BZN(l) 
normally equals 0) . 

Logical variable; if TRUE, the time correction 
for Briggs' formula is included in order to 
account for the instantaneous cloud formation. 
Ratio of ambient potential pressure to altitude. 
Time interval for cloud rise calculations. 

Note that this "DT" is different from the "DT" 
in the main program. 

Cooling effect of the deluge water on the bottom 
cloud . 

Entrainment constant. 

This logical variable if TRUE, causes additional 
printouts in the CLDRIS routine. 

Number of CLDRIS calculation intervals between 
CLDRIS printouts. 

Number of clouds. 

Ambient pressure, wind velocity, and temperature 
arrays, respectively. 

Heat content of the exhaust per second of fire. 
Time of fire for bottom cloud. 

Values used to calculate time of fire for upper 
clouds; t = TFA x (ALTITUDE) **TFB . 

Altitudes at which pressure, wind velocity, and 
temperature are given in P, U, and T arrays. 
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Ill . EXPLANATION OF IMPORTANT VARIABLES 


Name 


Symbol 


Description 


MAIN PROGRAM 


AN, BN, B 
ASLAB, ASLB 


DAN, DBN, DB 

4 nV B , 4H B , 

X X* 


AB. 

X 

DT 

At 

F 

F. 

1,3 


G 

G. 

X 

H 

h 

IP 


L 


N 


NCLT 


PMB 



RABSED, RBSED 


RAD, RM 

RL, RLP, RLPP, L . . , L ! . , L'.' . 

ij ij 13 

RMAA 

SIG L. 

TLIM, TOLD 1 


These are temporary storage arrays for the mass 
distributions. They correspond to ANP, BNP, and 
BNEW, but they only hold the values for one cloud 
These arrays keep track of the total number of 
particles lost through sedimentation for each 
mass class, over the entire cloud (material A/B 
and B, respectively). 

Change in number A/B, number B, fraction B in 
A/B for a given mass class, for a given cloud. 


Length of (main routine) time interval. 

This array is used to calculate how a particle 
formed in a collision is reclassified, i.e., 
if 2 particles, m. and m. 

13 

the new particle (mass = nu + m.) is split up: 


(m. > m.) collide, 

13 


F..*m. into mass class i; (1-F..)*m.., into 
13 x xj x+1 

mass class i + 1. 

Array of the mean thermal speed of the i'th 
size particle. 

Vertical thickness of the cloud. 

Pointer for printout, when N=IP, the mass dis- 
tributions are printed, and IP is incremented 
by IOINT. 


Counter for cloud number. 

Counter for time interval. 

When N=NCLT, DT is recalculated, and NCLT is 
incremented by NCLTI. 

Array storing total original mass of material 
B in each cloud. 

For a given mass class, rate of sedimentation 
loss (number /second/ml) for A/B and B, 
respectively. 

Arrays storing the radius and the mass, re- 
spectively, of each mass class particle. 
Collision rate arrays: RL, between A/B size i 
and A/B size j; RLP, between B. and B.; KLPP, 
between A/B^ and B j . 1 J 

Array storing total J original mass of material 
A in each cloud. 

"Mean free path" for particle values returned. 
Values returned from CLDRIS. TLIM is the time 
of the most recent CLDRIS calculation; TOLD is 
the time of the previous calculation. Values 
at time TIME are linearly interpreted between 
these two. 
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Name 

Symbol 

Description 

U 

y 

Viscosity of air. 

VLBEG, VLEND 


Cloud volume at beginning and end of (main 
routine) time interval. 

VSED 

V s.i 

Sedimentation velocity array. 

YK 

' v-j 
•rl 

Collision kernel array. 
CLDRIS ROUTINE 

AEMO 

m . 
air 

This array stores the mass of air in each cloud. 

AIRMO 

This is the amount of air in the 

cloud, from the previous calculation (from 

ARMO array) . 

AIRMP 

Mo 

Total mass of cloud at the beginning of the 
interval . 

AIRN 

m . 

air , new 

Mass of new air entrained in cloud. 

API, AVI 

P, V 

Arrays storing original pressure and volume, 
respectively, in each cloud. 

AZN 

Z 

c 

Stores the height of the center of each 
cloud . 

CPA 

Cp i 

Specific heat of each pollutant at 100° K 
intervals, from 200° to 1300° (set by DATA 
statement) . 

CPB, CPC 


Arrays calculated to provide easy calculation 
of the specific heat of a pollutant at a 
given temperature. 

CPD, CPE, CPF 


Calculated to provide the integral of the 
mass times the specific heat of a pollutant 
with respect to temperature, over any range 
from 298° to 1300° K. 

CPMIX 

C 

p,mix 

Specific heat of the cloud. 

FRACA 

F. 

1 

Rate of mass flow of each pollutant into 
the cloud from rocket exhaust (set by DATA 
statement) . 

I 

i 

Used in loops; points to one of the 6 
pollutants. 

IPR 


Pointer for output; when N = IPR, cloud 
parameters are printed, and IPR is incre- 
mented by IPINT. 

J 

j 

Refers to a temperature range of 100 s Kr: 
from 100* (J+l) to 100* (J+2). 

L 


Sub-cloud number. 

N 


(CLDRIS) time interval counter. 

PAMB, UAMB, 

^amb , ^amb , 
T 

Ambient pressure, wind velocity, and tem- 

TAMB 

perature, respectively, at the cloud 


amb 

altitude. 

PMAS 

M 

p 

Mass of pollutant I in cloud L. 

Q 

Q 

Heat content of each cloud. 
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QL 

BMW 

TFAR 


Name 


Symbol 


Description 


VOBS 

WN 


M 


wi 


V obs 

w 


PROR, DPR, VOL, 
DVOL, TEMP, 
DTK 


Energy loss due to adiabatic expansion. 

Molecular weight of each pollutant. 

Length of time interval over which rocket 
fires exhaust into each cloud. 

Volume of cloud. 

Increase in height of cloud in (CLDRIS) 
time interval. 

These arrays return values to the main 
routine for pressure, temperature, and 
volume in each cloud. PROR and DPR are for 
pressure, VOL and DVOL are for volume, 

TEMP and DTK are for temperature. PROR, 

VOL and TEMP contain values for the end of 
the interval, DPR, DVOL, and DTK return 
the change in the quantities over the 
interval. 
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IV, SAME! E PROBLEM INPUT DATA 


CASE 1 


CONSTANT PRESSURE, TEMPERATURE, E VOLUME (PTVC OPTION) TESTING 
S JIM 

PTVC*.TRUE.,CNV0L«l.O,TEND*l.E10,CNTEM*3O0.0,CNPRES*l.,NLAY*l, 

KKC*.TRUE.,XKC0N*6.E-10, 

CN0NX*1.E-20,BNEW* 30*0.0, 

RHO- l.,CNH»l.E 10 »ANP» 30*0.0, 

RMO*5.23E-1A,ALFA*10.0,IR»10,NT»5000,TI*.01,IOINT*500, 

NClTI*50,NZS*0, 

BNP*.180,19,71»1.897E3,1.308E5,8.54fc5,1.312EA, 30*0.0, 

SEND 



CASE 2 


VARSI AL2D3 DISTRIBUTION* DEBRIS -100*AL203 HASS 

$JIM 

BNP=1 .944E7 ,2. 11E7,8.90?E7,1. 388E7,7.4o4E5,4.022E4,2153. ,115.8, 
d. 23o, 15*0.0, 

AHP = 2»0.0,9.595E7 ) 9.354E7, 1» 13E7, 1.3o5E6, -1.648E5, 1.991 EH, 

2409. ,290. 3, 35.07,4.23o,. 51 17,. Ool 8,7. 466E-3, 10*0.0, 

NZS=0 , 

R M0=1. E- 17, ALFA=10.0,IR=1o,RH 0=3.0, N7=5000, 

7I=. 1, IOI NT = 1 00, TEND=525. 077, NCLTI=10,BNEW =30*0.0, 

$END 

$JIM1 

TF=17.0,TFA=.635,TFB=.484,DPHZ=.001098, 

BZN = 0 . 0 , o5 6 . 0 , 

D7= 1 . 0 , NLAY = 1 , IPIN m = 25 , 

ZP=0.0, 141 .4d,304.88,4?1 .04, 589.o3, 009.76, 847. 26, 9 14. 63, 970. 73, 1056.7, 
1219.51 , 1524.39, 1542.99, 1829.27,2051.82,2134. 15,2203.96,2331.4,2439.02 
Z7=0.0, 141 .46,304.88,471.04,589.03,609.76,847.26,914.03,970.73,1056,7, 
1219.51 , 1 524.39, 1542.99 , 1829.27,2051 .82,21 34. 1 5,2203. 96,2 33”' .4,2439.02 
ZU=0.0, 141 .4o,304.88,471 .04,589.63,009.76,847.26,914.63,970.73, 1056.7, 
1219.51 ,1524.39,1542.99, 1829.27,2051.82,2134. 15,2203.96,2331.4,2439.02 
P=1 .002,.987, .9o9, .951 , .938, .936, .91 , .903, .898, .888, .872, .841 , .839, .81 1 , 
.79, .782, .776, .764, .754, 

7=300.5, 298. 2, 29o.o, 295. 1 ,294.7,294.6,293.4,292.5,291.9,291 .2,289.6,287.2 
287. 1 ,285.3,284. 1,283.5,283.0,283.7,284.2, 

U=2. 57, 3*3. o,2*3. 08, 5*3. 6, 2*3. 08, 3. 6, 2*4. 63, 5. 15,o. 18,7.72, 

SEND 


V. SAMPLE OUTPUT 


CASE I 

CONSTANT PRESSURE* TEMPERATURE* £ VOLUME (PTVC OPTION) TESTING 

AEROSOL PROGRAM 

DENSITY ClF PARTICLE MATERIAL » 1.000 G/ML 

NO. OF PARTICLE CLASSES » 10 

TIME STEP SIZE INDEX - .0100 

MASS CF SMALLEST PAkTICLE • .52300E-13 GM. 

RATIO OF CONSECUTIVE MASS CLASSES (ALFA) - 10.0000 

NO. OF TIME STEPS * 5000 


LAYER » 1 


INITIAL 

MASS OF A/B « 

0. G. 


INITIAL 

MASS OF 8 • 

522.110 G 


INITIAL 

MASS ■ 522. 

110 



SEDIMENTATION 

VELOCITIES 


CLASS 

MASS* 

PG DIAM* MICRO-M 

VELOCITY. CM/SEC/ 

1 

• 52300000E 

-01 . 46396088 

.872355E-C3 

5 

523.00000 

9.9957392 

.302669 

10 

52300000. 

963.96088 

IB 7. 083 


SAMPLE COAGULATION CONSTANTS 



K ( I* J) • 

1.0E-10 ML/PARTICLE/SEC 

CLASS 

1 

5 10 

1 

6.000000 


5 

6.000000 

6.000000 6.000000 

10 

6.000000 

6.000000 6.000000 



to 


TIME • 0.0000 

TIME INTERVAL « 0 

DELTA TIME » 0.00000 

CLOUD 1 


PARTICLE SIZE and COMPOSITION SPECTRUM 


CLASS 

MASS* 

0 1 AM » 

NUiBI, 

FRACTI 


PG 

MICR0-M 

PER ML 

8 IN 

1 

. 52300E-01 

. 463 96 

0. 

0. 

2 

.52300 

.99957 

0. 

0. 

3 

5.2300 

2.15352 

0. 

0. 

A 

52.300 

4.03961 

0, 

0. 


5 

523.00 

9.95573 

0 . 

0 . 

6 

5230.0 

21.53516 

0 . 

0 . 

7 

52300. 

46.39609 

0 . 

0 . 

8 

. 523C0E + 06 

99.95734 

0 . 

0 . 

9 

• 52300E+07 

215.35156 

0 . 

0 . 

0 

. 52300E + 08 

463.96088 

0 . 

0 . 


TIME ■ 135072.6961 
TIME INTERVAL ■ 5000 

DELTA TIME ■ 117.63110 

CLOUD 1. 


o o o o 


**++*+*++*+* SEDIMENTATION ***♦♦**♦♦**♦ 

N ( B ) » LOSS (N0/ML/SEC) ACCUMULATED (NO) 

PER ML A/B 3 A/? B 

. 1800 
19.71 
1 E 9 7 . 

• 130 &L + 06 


.8540E+06 

.1312E+05 



PARTICLE SIZE AND COMPOSITION SPECTRUM 


************ sedimentation ************ 


CLASS 

MASS * 

DIAM, 

N(A, B)» 

FRACTION 

MB)* 

LOSS 

(NO/ML/SEC > 

ACCUMULATED (NO) 


PG 

MICRO-M 

PER ML 

B IN A/3 

PER ML 

A/B 

8 

A/B 

B 

1 

.52300E-01 

. 96396 

0 . 

0 . 

•1338E-03 

0 . 

•9059E-17 

0. 

•5099E-09 

2 

.52300 

.99957 

0 . 

0 . 

•1136t-01 

0 . 

.39956-19 

0. 

.222 2E-01 

3 

5.2300 

2.15352 

0 . 

0 . 

1.096 

0 . 

.1631E-11 

0. 

9.177 

A 

52.300 

9.63961 

0 . 

0 . 

85.92 

0 . 

.56756-09 

0. 

3002. 

5 

523.00 

9.99573 

0 . 

0 . 

9329. 

0 . 

. 1 309E— 06 

0. 

.2150E+06 

6 

5230.0 

21.53516 

0 . 

0 . 

•1399E+05 

0 . 

•1879E-05 

0. 

•5770E+06 

7 

52300. 

96.39609 

0 . 

0 . 

5876. 

0 . 

.35806-05 

0. 

.91906*06 

6 

. 52300E+06 

99.95739 

0 . 

0 . 

261.0 

0 . 

•6350E-06 

0. 

.91726*05 

9 

• 52300E+07 

215.35156 

0 . 

0 . 

1.031 

0 . 

•7902E-08 

0. 

399.7 

10 

. 52300E+08 

963.96068 

0 . 

0 . 

. 297 5E-03 

0 . 

,55676-11 

0. 

.1793 


NO. OF A/B PARTICLES - 0. NO. /ML 

NO. OF 8 PARTICLES ■ 29092.607 NO. /ML 

ALL MASSES IN MICRO-G/HL 
MASS OF A/8 PARTICLES ■ C. 

MASS OF b PARTICLES » 522. 06170 


MASS FRACTION OF B IN A/8 » I 

ACCUMULATED SEDIMENTATION LOSSES (G> 

MASS OF MATERIAL B LOST > . 98675G06t-0i (G> 

MASS OF MATERIAL A LOST > 0. 

TOTAL MASS LOST > .93675006E-01 

SEDIMENTATION LOSS RATE (G/SEC) 

A/B PARTICLES • 0. 

B PARTICLES • .570B3389E-06 

H» VOLUME AT BEGINNING ANO END OF INTERVAL • . 10000000E +11 1.0000000 i. 0000000 


VO 

U) 



VO 

■p- 


CASE 2 


VARSI AL 203 0 1 S TR IB UT I ON. DEBR IS *100* AL 203 MASS 
CRAM CALLED ' 

RhW ( I ) * 18.016000 

RMW ( I ) ■ 36.A65000 

RMW ( I ) * AA.D5A003 

RMWCIJ- 101.96100 
RMW ( I ) ■ 29.013A00 

RMW ( 1 1 • 159.69200 

UAM8» ASPS. RTS * 0. 

2.2623 2 • 598A01E-02 

A. TIMED* .522985A8E+13 .99510503 

TIMED- . 9898A565E-01 

INTERMEDIATE OUTPUTl NO I T> CM1 F , C Z.C 1> CP D» CP £ ( it Ml ) 

CPF(L»M1).T01.*T02.TNEW.TF1»TF2»F1»F2 

FNEWjRKP 

21 -.111136 + 16 • 3 5929E + 13 -.A3768E+08 -.13310E + 11 .A1897E+08 

9288.2 310. A3 310. A3 310. A3 

-• 19769E+13 .27193E+13 .27193E+13 -.19769E+13 

.672 A1E+07 .270A5E-05 

************************************************ 

TIME INTERVAL ■ 1 

LAYER » 1 

NEW TEMPERATURE ■ 310. A3 

TIME - 1.0000 

DELTA TIME - 1.0000 

MASS OF NEW AIR INGESTED INTO CLOUD ■ .50536E+11 

MASS OF OLD AIR IN CLOUD ■ 0. 

NO. OF ITERATIONS - 21 

OBSERVED VOLUME « .AA663E+08 

HEIGHT OF CLOUD • AA9.18 

CLOUD VELOCITY - A9.185 

GL ■ 0 • 

ADDITIONAL VARIABLES 2 ZM.C°MIX.GAM>Q(L).PAM8» T A M 6/ PM T. P I , VI , A I RMT ( L ) 

• 32 93 1 AE-01 0. 0. .1987555 + 12 .953368 

295.297 .182367E+09 .958696 .3153966+08 0. 

INTERMEDIATE OUTPUT « N0IT»CM1F.CZ.C1.CPD>CPE ( L . Ml ) 

CPF(L>M1)>T01»T02»TNEW>TF1.TF2.F1»F2 

FNEW. RK P 

21 -. 13289E+16 .A3255E+13 -.A3766E+08 -.13310E+11 .A1397E+08 

9288.2 308. 1A 308. 1A 303. 1A 

-.23321E+13 .6258AE+12 .6258AE+12 -.233215+13 

•62881E+36 .28209E-05 



2 


TIME INTERVAL • 

LAYER - 1 

NEW TEMPERATURE ■ 303.14 


TIME • 2.0000 

DELTA TIME « 1.0000 

MASS OF NEW AIR INGESTED INTO CLOUD » .10306E+11 

MASS OF OLD AIR IN CLOUD - .5Q536E+11 

NO. OF ITERATIONS - 21 

OBSERVED VOLUME ■ .558626+03 

HEIGHT OF CLOUD « 483.96 

CLOUD VELOCITY - 34.779 

QL ■ .1312Q7E+1G 

ADDITIONAL VARIABLES! ZM, C PMI X» GAM» 0 ( L ) > P AM B , T AM B. PMT» P 1 , VI. AIRMT { L ) 

• 344016E-01 .272567 1.33475 .1987556+12 .949593 

295.056 .1823676+09 .958696 .315396E+08 .149233E+14 

AEROSOL PROGRAM 

DENSITY OF PARTICLE MATERIAL - 3.030 G/ML 

NO. OF PARTICLE CLASSES ■ 16 

TIME STEP SIZE INDEX - .1000 

MASS OF SMALLEST PARTICLE - .10000E-16 GM. 

RATIO OF CONSECUTIVE MASS C L AS SE S { AL FA ) - 10.0000 

NO. GF TIME STEPS - 5000 


LAYER ■ 1 


INITIAL MASS OF A/B - .122936E+10 G. 

INITIAL MASS OF B * .121232E+08 G 

INITIAL MASS - .1241496+10 

SEDIMENTATION VELOCITIES 


CLASS 

MASS* PG 

DIAM, MICRO-M 

VELOCITY, CM/SEC/ 

1 

• 10000UOOE-04 

.13532772E-0I 

.402300E-04 

5 

. 10000000E+90 

.39927647 

.233856E-02 

10 

10000.000 

18,532772 

3.04753 

15 

. 10000000E +10 

860.21503 

699.450 


'O 

Ul 



CT\ 


SAMPLE COAGULATION CONSTANTS 



K(I,J) * 

1.0E-10 ML/PAR.TICLE/S EC 

CLASS 

1 

5 

10 

1 

15.32860 



5 

A 13 . 9299 

8.378233 

5A18.733 

10 

20796.19 

5A18.733 

6.16A520 

15 

3578971. 

• 5660900E+08 

.2616289E+10 


15 20 25 30 


• 5660900E +08 
• 2616289E+10 
6.1A215A 


************************************************************ 

TIME • 0.0000 

time interval ■ 0 


DELTA TIME » 0.00000 

CLOUD 1 


PARTICLE SIZE AND COMPOSITION SPECTRUM 


CLASS 

MASS, 

DIAM, 

N ( A » B ) , FRACTION N(B)» 


PG 

MICRO-M 

PER ML B 

IN 

A/B PER ML 

1 

.10000E-0A 

.01653 

0. 

0. 

. 19 A AE + 08 

2 

.10000E-03 

.03993 

0. 

0. 

•2110E+08 

3 

. 10000E-02 

.08602 

. 9595 E + 0 8 

0. 

. 89J7E+08 

A 

. 10000 £-01 

.18533 

.935AE+08 

0. 

•1388E+08 

5 

. lOOOOE+OO 

.39928 

. 11 30 E + 0 8 

0. 

.7A6AE+06 

6 

1.0000 

.86022 

. 1365E + 0 7 

0. 

•A022E+05 

7 

10. 000 

1.85323 

. 16A8 E + 0 6 

0. 

2153. 

8 

100.00 

3.99276 

, 19 91 E + 0 5 

0. 

115.8 

9 

1000.0 

8.60215 

2ACA, 

0. 

6.236 

10 

10000. 

18.53277 

290.3 

0. 

0. 

11 

1C0000 

39.92765 

35.07 

0. 

9. 

12 

• lOOOOE + 07 

86.02151 

A. 236 

0. 

0. 

13 

• 10000E+03 

185.32772 

.5117 

0. 

0. 

1A 

. 10000 E+ 09 

399.276A7 

.61 0 OE-9 1 

0. 

0. 

15 

.10000E+10 

860.21503 

.7A66E-32 

0. 

0. 

16 

• 10QG0E+11 

1853.27720 

0. 

0. 

0 . 

INTERMEDIATE OUTPUT* 

NOIT.CMIF 

•CZ,C1,CP0,CPE (L 

, Ml ) 


CPF C L 

»M1) ,T01,T02,TN 

EW»TF1» TF2, 

F 1, F2 



FNEw, 

RKP 





?1 

- • 15 72 1 E + 16 

. 5 1 A8 8 

E+13 -.A3768E+08 

-. 13310E+11 


9238.2 

306.11 

306.11 

306.11 

- 

• 2 725 1 E + 13 

• 7A636E + 12 

.7A636E+12 

"■ * 

27251E+13 


-.82083E+06 .26210E-05 


************ SEDIMENTATION ********* 
LOSS (NO/ML/SEC) ACCUMULATED (NO) 

A/B B A/B B 


• A1897E+08 



TIME INTERVAL ■ 500 

LAYER • 1 

NEW TEMPERATURE - 287.40 

TIME • 500.00 

DELTA TIME • 1.0000 

MASS OF NEW AIR INGESTED INTO CLOUD - .11988E+10 

MASS OF OLD AIR IN CLOUD « .45593E+13 

NG. OF ITERATIONS - 21 

OBSERVED VOLUME • .47230E+10 


HEIGHT OF CLOUD • 2124.1 

CLOUD VELOCITY • .20132 

QL ■ .293159E+11 

ADDITIONAL VARIA8LcS« ZM, C PM I X» GAM# Q ( L ) , P AM3 , TAMB»PMT»PI»VI»AIRMT(L) 
.344067E-01 ,24035b 1.39753 .198755E+12 .782974 

283.573 . 182367E+09 ,958696 .315396E+08 .130962E+16 


************************************************************ 
TIME ■ 522.2967 

TIME INTERVAL » 240 

DELTA TIME ■ 5.25077 

CLOUD 1 


vo 


I 



00 


PARTICLE SIZE AND COMPOSITION SPECTRUM 
CLASS MASS > DI AM* NU.B). FRACTION 




PG 

MICRO-M 

PER ML 8 

IN A/8 

1 


.10000E-0* 

.01853 

0. 

0. 

2 


.10000E-03 

.03993 

0. 

0. 

3 


. 10000E-02 

.08602 

• 2726E+05 

.1776 

* 


.lOOOOE-Oi 

.18533 

• 26 2*E+0 6 

.1*32 

5 


. 10000E+00 

.39928 

. 1008E+06 

.1069 

6 


1.0000 

.86022 

. 11 39E+0 5 

.*8696-01 

7 


10.000 

1.85328 

1099. 

.15976-01 

8 


100.00 

3.99276 

12G.1 

•6*85£-02 

9 


1003.0 

8.60215 

13.20 

. 3076E-02 

10 


10000. 

18.53277 

l.*0* 

.6185E-03 

11 


100000 

39.92765 

.1*8* 

.57876-03 

12 


. 10000 E+07 

86.02151 

.1572E-01 

.*53*E— 03 

13 


.100006+08 

185.32772 

.1552E-32 

. 3210E-03 

1* 


. 10000E+09 

399.276*7 

.12*86-33 

.19936-03 

15 


.10000E+10 

860.21508 

. 5968E-05 

•115*6-03 

16 


.100006+11 

1853.27720 

.3207E-07 

•9*19£— 0* 

NO. 

CF 

A/8 PARTICLES 

■ *03086.29 

NO./ML 


NO. 

OF 

B PARTICLES 

■ 28697.606 

NO./ML 



ALL MASSES IN MICRO-G/ML 
MASS OF A/B PARTICLES • .13921375 

MASS OF B PARTICLES ■ .18818338E-03 

MASS FRACTION OF B IN A/B • .16807E-01 

ACCUMULATED SEDIMENTATION LOSSES (G) 

MASS OF MATERIAL B LOST • 1*2091.33 (G) 

MASS OF MATERIAL A LOST ■ .58Q63631E+09 

TOTAL MASS LOST • .580778*06+09 


SEDIMENTATION LOSS RATE (G/SEC) 
a/s particles ■ 28***6.3i 

B PARTICLES > .51176011E-02 


************ SEDIMENTATION 


N(B)» 

LOSS (NO/ML/SEC) 

ACCUMULATED (NO) 

PER HL 

A/B 

B 

A/B 

B 

0. 

0. 

0. 

0. 

.1552E+12 

5.**1 

0. 

•2783E-08 

0. 

,202*E+13 

. 1620E+05 

•33A1E-0* 

•1986E-0* 

•26A5E+15 

•17*9E+15 

•1201E+05 

. 8667E-03 

•3967E-0* 

. 3603E+16 

.20666+15 

*82.0 

. 1071E-02 

•5117E-05 

•35905+16 

•278*E+1* 

3.695 

. *566E-03 

»l*82E-06 

. 1A9AE+16 

•1*70E+13 

•2020E-02 

. 1837E-03 

•337*E— 09 

•6291E+15 

•5758E+11 

. 6905E-08 

•8818E-0* 

• 5071E-1* 

• 309*E+15 

.38066+10 

0. 

.*3826-0* 

0. 

• 1570E+15 

.17376+09 

0. 

•2138E-0* 

0. 

.78796+1* 

130.9 

0. 

• 10*3E-0* 

0. 

•3986E+1* 

0. 

0. 

. 3815E-05 

0. 

• 1556E+1* 

0. 

0. 

.11106-05 

0. 

.52226+13 

0. 

0. 

.2083E-06 

0. 

.13**6+13 

0. 

0. 

.20*9E-07 

0. 

.26*26+12 

0. 

0. 

•2*68E-09 

0. 

• 88366+10 

0. 


H# VOLUME AT BEGINNING AND END OF INTERVAL ■ 


208399.*! 


*7390261E+10 


*73959086+10 


APPENDIX B 


PROGRAM LISTING 


1 


10 


15 


20 


25 


30 


PROGRAM CRAM( I NP JT, OUT PUT, PUNC H, TAPES, T AP E5- IN PUT, TAP£6»0UTPUT, 

1 TAPE7-PUNCH) 

DIMENSION AN(30),ANP(30,20), ASLAB( 20,30), ASLB( 20,30), BC 30), 

1 BN(30),BNEW(30,20),BNP(30,20),D(30),DPR(20),DTK(20),DVDL(20), 

2 F(30,30),G(30),PM(20),PMAA(20),PMB(20),PROR(20)«RAD(30), 

3 RL(30»3Q)»RLP( 30,30) »RLPP(30»30)»RM( 30 ),SIG13Q)» 

4 TEMPI 20), TITLE ( 79) , VOL ( 20 ) , VSED(30), XK< 30,30) 

DIMENSION FACTRA(20),FACTRB(20) 

LOGICAL PTVC, KKC 

COMMON /CA/CN3MX,NLAY, VOL, DVOL, TEMP, DTK, 

I PROR, DPR 

NAMELIST /JIM/ ALFA, ANP,BNEW, BNP, CNH, C NOMX, CNPRES, CNT£M»CNVQL, 

1 FACTRA, FACTR8, 

2 ICINT,IR,KKC,NCLTI,NLAY,NT,NT1,NT2,NZS,PTVC,RM0,RN2M, 

3 PHD, TEND, TI,XKCDN 

C THIS PROGRAM IS INTENTED TO CALCULATE THE COAGULATION OF TWO TYPES OF 
C PARTICLES IN THE EXHAUST CLOUD OF A ROCKET BLAST OFF. THt PARTICLES ARE 
C DIVIDED INTO A NUMBER OF DISCREET MASS SIZES. 

C JIM IS THE INPUT DATA NAMELIST FOR THE MAIN PROGRAM } THE VARIABLES ARE« 

C RMO « MASS OF THE SMALLEST PARTICLF (SMALLEST MASS SIZE) 

C ALFA ■ MASS RATID OF CONSECUTIVE MASS SIZES 
C IR ■ NUMBER OF MASS SIZES 
C RHO ■: DENSITY OF PARTCILt MATERIAL 

C NT ■ NUMBER OF TIME INTERVALS TO BE CALCULATED, IF CLOUD RISE TIME (TEND) 

C IS NOT EXCEEDED FIRST 

C ANP,3NP,bNEw • ARRAYS GIVING A/3 DISTRIBUTIONS, B DISTRIBUTIONS, AND 
C FRACTION a IN A/B, RESPECTIVELY. EACH ARRAY HAS 600 ELEMENTS J 1-30 FOR 

C THE BOTTOM CLOUD, 31-60 FOR THE 2'ND CLOUD, UP TO 20 CLOUDS; ELEMENTS 1, 

C 31, 61, ETC. GIVE NUMBER FOR SMALLEST MASS PARTICLE IN EACH CLOUD. 

C TI » INITIAL VALUE FOR TIME STEP SIZE INDEX, 

C IOINT - NUMBER OF TIME INTERVALS BETWEEN PRINT OUTS. 

C TEND » TIME 4T WHICH TO STOP CALCULATIONS (CLOUD RISE TIME), IF NT NOT 
C EXCEEDED FIRST. 

C NT 2 « NUMBER OF TIME INTERVAL AT WHICH TO FIRST REVISE TlJ THEN COUNTER 


VO 

VO 



100 


35 


40 


45 


50 


55 


60 


65 


70 


C FOR SUBSEQUENT REVISIONS OF TI. 

C NT1 - INCREMENT VALUE FOR NT2 WHEN TIME INTERVAL NT2 REACHED 
C RN2M • REVISE VALUE FOR TI < TI*TI*RN2M> 

C NCLTI • NUMBER OF TIME INTERVALS INBETwEEN RECALCULATIONS OF DELTA TINE 
C <'DT'). 

C END OF NAMELIST VARIABLES. 

C FOR ANY DO LOOPS USING N ,L* OR I* N IS THE COUNTER FOR THE TIME INTERVAL* 

C L POINTS TO THE CLOUD* AND I IS THE MASS SUE POINTER. 

C IP IS THF COUNTER FOR THE PRINTOUT* IE* WHEN IP»N* MASS DISTRIBUTIONS ARE 

C PRINTED* AND WHEN THE TOP LAYER HAS BEEN PRINTED* IP IS INCREMENTED BY 

C 'IOlNT'. 

C 'NCLT' IS THE COUNTER FOR RECALCULATIONS OF *DT'* WHEN NCLT*N, *DT' IS 
C RECALCULATED* AND 'NCLT' IS INCREMENTED BY 'NCLTI*. 

C 'TIME' POINTS TO THE BEGINNING OF THE INTERVAL BEING EXAMINED. 

C ' CLDRIS ' IS A SUBROUTINE THAT CALCULATES TEMPERATURE* PRESSURE* AND 
C VOLUMEOF EACH CLOUD. THE FIRST TIME IT IS CALLED AS 'CLDRIS'* AFTER 
C THAT* AS 'POINT*. WHICH IS AN ENTRY POINT IN 'CLDRIS'. 'CLDRIS* SETS 
C VALUES IN 6 ARRAYS* • VOL* TEMP* PROR* DVOL* DTK* DPR • . THE FIRST THREE GIVE 
C THE VALUES OF VOLUME* TEMPERATURE* AND PRESSURE* RESPECTIVELY* IN EACH 
C CLOUD* AT TIME ' TLIM* • THE SECOND THREE GIVE THE CHANGES IN VOLUME* 

C TEMPERATURE, AND PRESSURE, RESPECTIVELY, IN EACH CLOUD, SINCE THE LAST 
C TIME THEY WERE CALCULATED ( 'TOLD* ) • ALL THREE CHANGE ARRAYS SHOULD CONTAIN 
C POSITIVE VALUES* ASSUMING DROPS IN TEMPERATURE AND PRESSURE* AND A RISE 
C IN VOLUME. WHENEVER 'TLIM' IS EXCEEDS BY 'TIME'* 'POINT* IS CALLED, 

1 READ (5*111 ) (TITLEm*I»l*79) 

DO 2 1*1*20 

FACTRAt I)*1.0 
F ACTRB < I ) *1 .0 
Du 2 J*l» 30 
ASL AB ( I* J ) *0,0 

2 ASLB ( I * J ) *0.0 
111 FORMATtlX* 79A1 1 

Wfc ITE (6*111) (TITLE(I)*I*1*79) 

NT2-1000 

N2S»1 

NT1-1000 

RN2M*1.0 

NLAY-1 

PTVC*. FALSE. 

KKC*. FALSE. 


TOT 


75 


30 


85 


90 


95 


100 


105 


110 


READ <5, JIM) 

00 112 J«1*NLAY 
00 112 I-l * IR 

ANP(I,J>«ANP( I* J)*FACTRA(J) 

112 8NP( I* J)»8NP( I* J)*FACTRS( J) 

IP * 10 INT 

TKOLD-IOOOOO.O 

T1PER-TEN0/103.G 

NCLT-1 

TIME-0.0 

C FOR ? TVC > WHEN THIS IS SET TRUE BY THE INPUT DATA (DEFAULT VALUE- 
C FALSE )» CONSTANT PRESSURE* TEMPERATURE* AND VOLUME ARE ASSUMED. GIVEN 
C IN INPUT AS CNPRES ( ATM. ), CNTEM(K)* CNV0L(M*+3)» ALSO* NLAY MUST BE 
C GIVEN IN THE INPUT OATA BECAUSE CLDRIS IS NOT CALLED. CNOMX MUST BE 
C GIVEN, ANO HEIGHT(CNH) MUST ALSO BE GIVEN. (HEIGHT IN CM) 

IF(PTVC) GO TO 80 
CALL CLDRIS(TLIM) 

SO TO 81 

80 TLIM-TEND+1.0 
DO 82 I-l.NLAY 
VOL ( I ) -CN VOL 
DVOL ( I ) -U. J 
PROR ( I ) -CNP RE S 
DPR (I )-U.O 
TE MP ( I ) -CNTEM 
32 DTK(I)-0.0 

bl TOLD-O.O 

DO 30 N * 1 * N T 
1F(N.NE.NT2) GO TO 11 
NT2-NT2+NTI 
TI-TIWRN2M 

C FOLLOWING IS THE CLOUD LOOP* IE* EVERY TIME THROUGH THIS LOOP* THE 
C DISTRIBUTIONS FOR ANOTHER CLOUD ARE CALCULATED. 

C ' TK • AND *P* ARE VALUES OF TEMPERATURE AND PRESSURE* INTERPOLATED FOR 
C TIME 'TIME'. *U* IS THE VISCOCITY OF AIR* • RLAMDA * IS THE MOLECULAR MEAN FREE 

C FREE PATH. *ALFC* IS A DUMMY VARIABLE EQUAL TO THE MASS OF THE PARTICLE 
C IN QUESTION. *RK* IS BOLTZMANN* S CONSTANT. 

11 DO 60 L-1*NLAY 


ii5 


120 


125 


130 


135 


1*0 


1*5 


130 


TK-TEMP (L)+ (TLIM-TIME ) / (TL IM-TOL D ) *DTK < L ) 

P«PROR(L)+(TLIK-TIME )*DPR< U / (TL IM-TOL D) 

U-1.055E-5+S0RT (TK) 

RL AMDA-2. 2E-B+TK/P 
ALFC*1. /ALP A+RMO 
RK*1. 36E-16 
Vi»RK*TK/lu.85/U 
V2-1.25 7+RL AMDA 
V3«. **»LAMQA 
V* — i.I/RLAMDA 

C FOLLOWING IS A LOOP TO CALCULATE VALUES DEPENDANT ON MASS SIZE. ' AN# BN* AND 
C B' ARE USED TO STORE THE MASS DISTRIBUTIONS FOR THE LAYER BEING 
C CALCULATED# ALLOWING NcW VALUES TO BE STORED IN *ANP»BNP# AND BNEW 1 AND 
C STILL HAVE THE OLD VALUES AVAILABLE. 

C • R M ( I ) • « THE MASS# 'RAD(I)' • THE RADIUS# AND 'VSEO(I)' • THE 
C SEDIMENTATION VELOCITY OF A PARTICLE OF MASS CLASS I. 

C * S I G ( I ) • ■ VALUES USED TO CALCULATE THE COLLISION KERNEL ARRAYS. 

DO 10 I*1#IP 
AN ( I ) • ANP ( I , L ) 

BN ( I )-BNP(I,L) 

3(I)*SNEW(I#L) 

AL FC-AL FC + ALF A 
RM ( I ) -ALFC 

RADI- (.238 7*ALPC/RH0)**(1. /3. ) 

RADII ) - RAD I 
TE S T- V* *R AD I 

IF (TEST. LE. -500. ) TEST— 500.0 

01 -VI /RAD I M i.+V2/RADI+V3/RADI*EXP{TEST) ) 

D( I ) • D I 

GW.5A0+RK + TK/ ALFC 
Oil )-GS 
GR-SQRT(GS) 

DUM-2.5*6+DI/GR 

SIGH )«.09256*GR/DI/RADI*{ (2«*RADI+DUM)**3.0-(*.*RADI+RADI 
1 +DUM+DUi)**1.5)-2.326*RADI 

10 IFIRADI.LE.2.0E-3) VS ED 1 1 ) -980. 0* AL FC + D I / RK/TK 

IF(.96.LE.TK/TK0LD.AND.1.0*.GT.TK/TK0LD) GO TO *3 

TKOLD-TK 

DO ** I»1»IR 

RADI -RADI I ) 

IFIRADI.LE.2.0E-3) GO TO ** 



155 


160 


165 


170 


175 


130 


185 


190 


TERM1-2.74EZM (P*RAD I ) *♦0.666667 )/TK 

TtRrt2»4.93£fc*RM(I)/(SaRT(TK) )/RADI 

EN-TERM2/2D.0 

EP-TERM2 

DO 47 J*l,25 

VS£DI«{£N+cP)/2,G 

FNE W» VSED I+TE RM1*(VS EDI **1. 666667 )-T£RM2 
IFLFNEW.GT.O.O) GO TO 48 
EN-VSEDI 
GO TO 47 

43 EP-VSEDI 

47 CONTINUE 

VSED(I)»VSEDI 

44 CONTINUE 

43 4«0.0 

C The FOLLOWING LOOP CALCULATES THE F ARRAY* THE *XK» (COLISSION KERNEL) 

C ARRAY, ANO SOME OF THE VALUES USED FOR CALCULATING ‘DT*. 

00 20 I ■ 1» I R 
00 20 J-1,1 
RSUM«RAD(I)+RAD(J) 

Q$UM»D( I)+D ( J ) 

0UH«1.+RAD< J)/RAD(I) 

X< I * ( DUM*0JH-1 . /DUN) *3«1416*PAD{I)*RAD(I)*(VSED(I)-VSED(J))+ 

1 12.57*RSUM*DSUM/<RSUH/(RSUM*SQRT(SIG(I )*SI3(I)*SIG( J)*SIG(J) )) 

2 +4*DSUN/ RSUM /SORT (G(I)+G(J))) 

IF(KKC)XKI»XKCON 

<K(I*J)»XKI 
XK ( J » I ) "XK I 

F ( I, J)«l-R1(J)/( ALFA-1.0) /RM( I) 

IF (L.Nt.l.QR.N.NE.NClT) GO TO 20 
ATT»XKI*RN( J) /RM( I) 

A«AHAX1(A*ATT*AN( J), ATT*BN( J) ) 

20 CONTINUE 

C THE FOLLOWING (UP TO STMT 28) CALCULATE SOME OF THE INITIAL MASSES* 

C AND PRINT OUT THE INITIAL CONDITIONS, »H* ( THE THICKNESS OF THE CLOUD) 
C IS SET AT lQJ METERS* EXCEPT FOR CLOUD 1* WHERE IT IS SET AS THE RADIUS. 
C »DT* IS CALCULATED, IF »N«NCLT», AND *L»1'. 

IF < N.Nt • 1 ) GO TO 28 
PM { L ) «G»0 



104 


! 


195 


200 


205 


Z10 


215 


220 


225 


230 


PMB8*0.0 
P^A3*C. 0 
PMAA1L)*0.0 
DO 25 I - 1 * IR 
PMBB*PMB6+RM< I )*3N(I ) 

PflAB=PMAB+RM(I)*AN(I ) 

PflAA!L)»PMAA(L)+RK(I »*AN(I )*(1-B(I) ) 

2 5 FM(L)*PM(L)+RM(I)*(AN(1)+8N(I)) 

PMB(L)*PM1L)-PMAA(L) 

H100-H/100 

IF (L.EQ. 1) WRITE 16. AOO)RHO» IR»TI,RMQ,ALFA,NT 
WRITE ( fa* A05 ) L 
A05 FORMAT!* LAYER « *15) 

DUMMVOt (L)-DV0L(L))*1.0E6 

PMAA(L>*PMAA!L)*OUM 

PMB(L)*PMB(L)+DUM 

WRITE <6, AOo) PMAB+OUM, PMB8*DUM,PM(L)*DUM 
A06 FORMAT!//* INITIAL MASS OF A/B ■ *G1A.6* G.*/ 

1 * INITIAL MASS OF B « +G1A.6* G*/ + INITIAL MASS ■ +G1A.6) 
WRITE (6» A 10) 

1*1 

WRITE (6, A20)I,RM(l)+i.E12,RAD(l)+2.EA, VStD(l) 

DO 500 I»5#IR»5 
RMI*RM( I )*1.0E12 
kADI»RAD(I)*2.0EA 

500 WRITE! b»A2Q) I » RMI, R ADI » VS EO ! I ) 

wR ITE ( 6* A30 ) 

1-1 

WSITE(6>AA0) I,XK(1» 1)+1.E10 
DO 5iC I*5*IR»5 

510 W P I T E ( & » A AO ) I>XK(I>1)*1.E10»(XK(I>J)*1.E1Q>J*5»IR»5) 

AOO F0RMAT(20X*AER0SJL program*// 

2 * DENSITY OF PARTICLE MATERIAL « *F7.3* G/ML*/ 

3 * NO. OF PARTICLE CLASSES - +15/* TIME STEP SIZE INDEX ■ + 

A F9.A/+ MASS OF SMALLEST PARTICLE * +G12.5* GM.*/ 

5 + RATIO OF CONSECUTIVE MASS CL ASS ES ( ALF A ) » *F8.A/ 

6 * NC. OF TIME STEPS *+I5//) 

A10 FGPMAT!10X,*SEDI(1ENTATIQN VELOCITIES*/* CLASS*I2X*MASS» PG* 

1 7X*DIAM, MICR0-M+AX*VE10CITY» CM/SEC /+ ) 

A20 FORMAT! IX, I5» 2 ( 5X» G15. 8 ) , 5X* G15. 6) 

A30 FGRMAT(//10X,*SAMPLE COAGULATION CONST ANTS + /15X* 
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235 


240 


245 


250 


255 


260 


265 


27 0 


1 33HK(I,J) * 1.0E-10 NL/PARTICLE/SEC/1X.*CLASS+14X*1*14X*5* 

2 13X410413X*15*13X*20*13X*25 + 13X*30*/ ) 

440 F0RMAT{lX,I5,4Xf 7(1X»GA4.7)) 

NZERO-O 

*RITe(6,4'50) TINEjNZekO»DT>L 
DO 730 1>1j IR 

730 HR I T£ ( 6» 460 ) I , Rfi ( I ) 41 . OE 1 2, R AO ( I ) *2. 0E4» AN ( I ) , B ( I ) > BN C 1 ) 

2a H-100C0.0 

IHL.NE.l) GO TO 29 

H- ( ( VOL (L)-(TL IN-TIME )*DVOL(L)/< TLIM-TQLD) )+*(l. 0/3.0) ) *124. 07 

IF(N.NE.NCLT) GO TO 29 

NCLT-NCLT+NCLTI 

OT-TI/A/IR 

IFtDT.GT.TIPER JQT-TIPER 

IF (DT.GT.T1PER/3.0I IOINT«INT < T1PER*5 /DT > 

C THE COLLISION RATE ARkAYS ARE CALCULATED. 

29 DO 50 I«1»IR 

DO 50 J-I.I 
XKI«XK(I,J ) 

DUh»XKI*AN(I)*AN(J) 

RL( I » J ) "QUM 
RL( J> I )"DUH 
DUH>XftI*BNm*BN( J) 

RLP ( I. J ) *DUM 
R LP < J. I ) *DUM 

RLPPCI»J ) »XKI *AN ( I )*BN l J ) 

50 RL PP ( J* I )»XKI*Ai'((J)*BN(I) 

C THIS IS THE MAIN LOOP FJR CALCULATING NEW NASS DISTRIBUTIONS. 
VLEND«VOL(L)-DVOL(L) +( TLIM-TIME-DT ) / ( T LIN-TOLD ) 
VLBEG"VCL(L)-DVOL(L)*<TLIN-TINE)/(TLIN-TOLD) 

IF (PTVC ) H-CNH 
DO 100 I»1»IR 
D AN»0. 3 
DBN"0.0 
OB»C.O 
IM«I-1 

IF(I.EQ.l) GO TO 55 
PMM-RMIH) 

DUM«B<IH) 

RMNbN"RNN*DUrt 
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275 


280 


285 


290 


295 


300 


305 


310 


55 P MI *RM ( I ) 

RM6I«RMI+B <I» 

C THIS LOOP PERFORMS THE SUMMATIONS FOR EACH MASS SIZE 'I'. 
DO 200 J«1,IR 
S-l. 0 

IF(I.EQ.J) S».5 
SM-1.0 

IF (IM.EQ.J) SM»0.5 
RLV-RL(I> J )+DT 
PLPV*RLP(I, J)+DT 
RLPPV«RLPP(I> J)*DT 
P LP PR*RL PP ( J» 1 >*DT 
DBN»P6N-RLPV-RIPPR 
DUM-RLV+RLPPV 
DAN-DAN-OUM 
DB*DB-DUM+RMB1 
IF (J.GT.I) GO TO 200 
PMJ«RM( J) 

RMBJ»RMJ+B<J) 

DUM-PMI/ (RMJ+RMI ) 

FV*F ( I* J ) 

SFV«S*FV 

SLFV«SFV*RLV 

FLPPV-FV+RLPPV 

DAN-DAN+SLFV+FLPPV 

DB«DB+SLFV*DUM*(RM9I+RMBJ)+FLPPV*DUM+< RMBI + RMJ ) 
DBN-DBN+SFV+RLPV 
IF(J.GE.I.OR.I.cQ.l) GO TO 200 
DUM2*RMI/(RMM+RMJ) 

F M V K 1-F ( IM» J ) 

SFMV«FMV*SM 

SIFMV«SFMV + RLUM> J)*DT 
FLPPM«RLPP<IM, J)*FMV+OT 
FLPPk*RLPPR*F V 
iiAN-DAN + SLFMV+FLPPM+FLPPR 

D8»DB + SLFMV*0UM2MRMMBM+RMBJ ) +FL PPM*0UM2+ ( RMMBM+RM J > 

1 *FLPPR+DUM+(RMI+RMbJ) 

DBN»DBN+SFMV+RLP(IM> J)*OT 

IF( J.GE.IM.0R.I.LE.2) GO TO 200 

FLPPMP»RLPP<J>IM)*FMV+DT 

DAN-DAN+FLPPMR 
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315 


320 


325 


330 


335 


340 


345 


350 


DB»UB+FLPPMR*DUM2*(RMN+RMBJ) 

200 CONTINUE 

C STORE NEW NASS DISTRIBUTIONS AN • ANP, BNP, BNew TAKING INTO ACCOUNT 
C THE INCREASE IN VOLUME, THE SEDIMENTATION LOSS, AND THE SEDIMENTATION 
C GAIN, (WHERE APPLICABLE). 

DUM«AN(I ) 

DUM2-DUM+DAN 
PUM3-BNU J+D3N 
IF(DUM2.LT.CNUHX) DU(12»0.0 
IF(DUM 3.LT.CNQMX) DUM3-0.0 

35 IF ( DUM2 • E a .0 . 0 ) GO TO 92 

P NE W ( I,L)»(RMI*DUM+3<1)+DB)/RMI/DUM2 
GO TO 95 

92 BNE W ( I, L ) “0.0 

95 VSEDIH»VSED(I )/H 

VSEGPR*VSEDIH*DT 
IF ( VSEDFR. LE. 1.0) GO TO 37 
WRITE (6,36)VSEDFR,DT,H,VSED< I) 

VSEDIH-l.O/DT 

36 FOPHaTI* SEDIMENTATION FRACTION .GT. It -+G14.6* DT» *G14.6 

1 * H« * G1 4 . 6+ VS ED ( I ) » * G14.6) 

37 RABSfcD«DUM2*VSEDIH 
ABSED-RABSED+DT 
R3SED-DUM34VSEDIH 
bSED»R3SED*DT 

DTL-(TLIN-TIHE-DT)/(TLIN-T0LD) 

ASLAB(L,I)«ASLAB(L»I J+ABSE0*VLBEG*1.0Eb 

ASLB(L,I)*ASL8(L,I)+BS£D*VLBEG*1.0E6 

IF(L.cO.l) GO TO 62 

DVL-VLEND/ ( VOL (L-l )-OTL*DVOL ( L-l ) ) 

ANP(I,L-l)»ANP(I,L-l)+ABSED*nVL 
BNP(I,L-1 J«BNP(I,L-1)+BSE0*DVL 
62 G ( I ) *R ABSED 
D(I)-PBSED 
DVL"VL3EG/VLEND 
ANP ( I , L ) * ( DUM2-ABSED ) *D VL 
BNP(I,L)«(DUM3-BSED)*DVL 
100 CONTINUE 

C PRINT OUT IF 1 N» N P • , OTHERWISE GO TO NEXT CLOUD. 

IF(N.Nt.IP) GO TJ 60 
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35 5 


360 


365 


370 


375 


360 


365 


390 


IF(L»EO.NLAY)IP"IP*IdlNT 

SLPAB-0.0 

5LRB-0.0 

ABNO-O.O 

BNO-O.O 

S MAB*0» 0 

SHB-o.0 

FNBINA-O.O 

VRITE(6»A50)TINE>N»0T*L 
7o0 DO 70 I • 1# I R 
RMI*RM< I ) 

SLRAB«SLkAB+R«1*G(II 

SLRB»SLR8+RMI*D<I) 

B I a BNE u ( I * L ) 

AN1»ANP(I,L) 

BNI"BNP ( I» L ) 

A3ND* ABNO+ ANI 

DUM" AM*RMI 

SMAB»SNAB+DUM 

FMBINA«FMBINA+DUM*BI 

BNO-BNO+BNI 

SMB a SMB+6NI*RNI 

RMI«RMI*1.E12 

RADI«RAO(I >*2.0EA 

WRITE(6>A60)I»Rf1I>RADI>ANI>BI»BNI>G(I)»D(II»ASLAB(L»I)#ASLB(L»I) 
70 CONTINUE 

SMBINA-FMBINA/SNAB 

BSEDL-PMB (L >-(FNBINA + SMB)*VLEND*1.0E6 

AS£DL»PMAA(L)-(SHAti-FMBINA )*VLEND*1.0E6 

SlRAB«SLRAB*1.0E6+VLEND 

SLRB»SLRB*1»0E6*VLEND 

SMB«SMB*1.0E6 

SMAB»SMAB*1.0E6 

TSEDL“BSEDL+ASEDL 

WRITE (6» 470 > ABN0.»BN0»SMAB»$MB»SM8INA* BSEDL» A$EDL*TSEOL» 

1 SLRAB»SLRB«H«VLb£G*VLEND 

*50 FGRrtAT(//IX#15(*H****>/* TINE ■ *F11.*/* TIME INTERVAL - * 

1 15/* DELTA TIME ■ * F12.5/* CLOUD* 

1 I5///10X,*PARTICLE SIZE AND COMPOSITION SPECTRUM*/ 

7 71X» 12 ( 1H* )» * SEDIMENTATION *#12(1H*)/* CLASS*6X*MASS**9X 
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395 


8 +DIAM»*5X*N(A»8)»*3X + FR ACTI0N*6X*N( B)**6X*L0SS (NO/KL/SEC)* 

9 7X * AC C U*UL AT ED t MU ) * / 1AX*PG*9X*MICRQ-M*5X*PER ML*AX*8 IN A/B* 

1 A X * P E P ML+»2(8X*A/8*10X*B*)) 

A6C FORMAT l IX, I3,5X,G12.5,F11.5,AX,7G11.A> 

A 70 FORMAT ( /* NO. OF A/9 PARTICLES » *616.6+ NO. /ML*/ 

AOO 1 • NO. OF B PARTICLES ■ +G16.8* NO. /ML*//* ALL MASSES IN * 

2 *M ICR 0-G/ ML */ + MASS OF A/B PARTICLES » *616.8/ 

3 * MASS OF B PARTICLES ■ +G16.8/+ MASS FRACTION QF B IN A/B • + 
A G12.5//+ ACCUMULATED SEDIMENTATION LOSSES (G) */ 

5 + MASS OF MATERIAL B LOST « *G16.8* (G)*/ 

AOS 6 + MASS OF MATERIAL A LOST • *616.8/* TOTAL MASS LOST • + 

8 G16.8//* SEDIMENTATION LOSS RATE (G/SEC)* 

9 /+ A/B PARTICLES » * G 16. 8/ * B PARTICLES » * 

7 G16.8//* Hj VOLUME AT BEGINNING AND END OF INTERVAL « * 

8 3G16.8///) 


A10 

60 

CONTINUE 

TIME-TIME+OT 

IF (TIME. GT. TEND) GU 

TO 

31 



IF (TIKc.LE.TLIM) GO 

TO 

30 


32 

TQLD-TLIM 



A16 

30 

CALL POINT (TL IM ) 

IF ( TIME.GT.TLIM) GO 
CONTINUE 

TO 

32 


31 

IF(NZS.NE.l) GO TO l 
STOP 



A20 


END 
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SUBROUTINE CLDRIS 


1 


5 


10 


15 


20 


25 


30 


3 5 


SUBROUTINE C LDR I S ( TI Mt ) 

DIMENSION CPA16, L2>, CPB< 6* 11), CPC(6*11)*CPD( 20* 11)*CP£(20*11>» 

1 CPF (20, 11 )* FRACA (6>,P(100),PMAS(20,6)*RMW(6),T1100)*U1100)* 

2 ZP(100),ZT(100)»ZUC100) 

DIMENSION AIRHT(20),API(20),APMT(20)*ARMG<2G>*ASPS(20)* 

1 AVI(2G)*AZN(2Q)*3ZN(20),DPR(20)*DTK(20)»DVOL(20)*PROR(20)* 

2 C(2G),RTS(20),TEMP(20),TFAR(20),VOL(20) 

LOGICAL CDRT-WINTPR 

COMMON /CA/CNOMX»NLAY*VOL*DVOL»TEMP*DTK* 

1 PRUR*DPR/Cd/ZT»T 
COMMON /CC/ZP* P/lG/ZU*U 

DATA FRACA/ 1. 39 7o E6* 9.Q168E5* 1.7551E6* 1.2626E6* 5. 0773 E6* 1. 665E5/, 

1 PMw/18. 016* 36. 565* 55. 055* 101. 961* 28. 0135* 159.692/ 

DATA ((CPA(I*J)*J=l»i2l*I*l»6)/ 

1 7.969*8.027* c.i. 2 0*8.515*3. 676*8. 955*9.256*9. 557*9. 851* 

2 10.152*10.555.10.723* 

3 6. 961*6. 9 65* 6. 973*7. 005, 7. 069, 7. 167* 7. 239, 7. 5 23* 7. 559* 

5 7.693,7.815* / • 9 36 * 

5 7.715* 8.3 9 6* 9.8 77,10.666*11.31, 11.656*12.293*12.667*12.988* 

6 13.253*13.566*13.656* 

7 12. 22 *13. 98, *22. 965 *2 5. 366, 26. 8 99* 2 7. 956* 28. 7 13, 29. 3 17* 29. 8 21* 

3 30.26*30.o53*31.008* 

9 6. 957, 6. 961, 6. 99* 7. 065, 7. 196* 7. 35* 7. 5 12. 7. 6 7* 7. 815* 

1 7.555*9.061,8.162* 

1 16. 3, 25. 9, 28. 71. 31. 5*33. 75*35. 786* 37. 815* 39. 792*36.0* 

2 33.656*33.322,33.998/ 

NAMELIST /Jill/ BZN*DPHZ*DT*DWATER*GAMTWQ*IPINT*NLAY*P*QC»TF* 

1 TFA*TF8»T*U*:p*ZT*ZU»C0RTM,INTPR 

C THIS SUBROUTINE CALCULATES THE TEM PER ATUR E* PRE S SUR E» VOLUME* AND HEIGHT 
C OF THE CLOUD AS IT RISES. THE TEMPERATURE, PRESSURE* AND VOLUME ARE 
C RETURNED TO THE MaIN PROGRAM. 

C J I Ml IS THE INPUT DATA NAMELIST; THE VARIABLES ARE: 

C TF <= TIME CF FIRE FjR BOTTOM CLOUD. 

C TFA * X VALUES USED TO CALCULATE TIME OF FIRE FOR UPPER CLOUDS; 

C TF3 - X TIME » T F A ♦ (ALTITUDE) EXP ( T F B ) . 

C BZN » ARRAY GIVING INITIAL ALTITUDES OF BOTTOM OF EACH CLOUD* AND TOP OF 
C TOP CLOUD (EZN(l) NORMALLY EQUALS 0) 

C D PHZ - RATIO OF AMBiANT POTENTIAL PRESSURE fO ALTITUDE. 
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C DT - TIME INTERVAL FOR CLOUD RISE CALCULATIONS. NOTE THAT THIS 'DT' IS 
C DIFFERENT FROM THE 'DT' IN THE MAIN PROGRAM. 

C NL AY - NUMBER OF CLOUDS. 

C ZP « / ALTITUDES AT WHICH PRESSURE. WIND VELOCITY. AND TEMPERATURE ARE 
C ZU » / GIVEN IN P, U. AND T ARRAYS RESPECTIVELY. 

C ZT - / 

C P - + 

C U - + AK3IANT PRESSURE. WIND VELOCITY. AND TEMPERATURE ARRAYS. RESPECTIVELY 
C T * + 

C END uF NAMELIST 

C THE FIRST PART OF THE PROGRAM. UP TO ‘ENTRY POINT*. IS EXECUTED ONLY 
C ONCE . 

C FOR ALL LOOPS IN THIS ROUTINE. *L‘ POINTS TO A CLOUD. *1* POINTS TO ONE 
C OF 6 POLLUTANTS. AND * J * REFERS TO A TEMPERATURE RANGE OF IDO DEGREES K. 

C FROM 10CJ+100. TO 100J +200 K. 

C QC IS ThE HEAT OF THE EXHAUST RELEASED PER SECOND. 

C THE FOLLOWING ONE DIMENSIONAL ARRAYS ARE INITIALIZED. DEPENDANT ON »L‘. 

C ' Q ' IS THE HEAT CONTENT OF EACH CLOUD. 

C ‘RTS* IS THE SQUARE ROOT UF ‘S‘» USING TEMPERATURE AT THE GROUND. FOR 
C BOTTOM CLOUD. OR AT THE BOTTOM OF THE CLOUD OTHERWISE. 

C ‘ASPS' IS A PARTIAL CALCULATION OF THE Z EQUATION. INCLUDING 'A'. 

C TEMP. AND PRESSURE (FOR RHO OF AIR) TAKEN AT GROUND OR BOTTOM OF CLOUD. 

C WIND VELOCITY ( L NOT EQUAL 1) AT MIDDLE OF LAYER. 

C 'TFAR' IS THE TIME OF FIRE FOR THE CLOUD. 

C 'AVI' IS THE ORIGIONAL CLOUD VOLUME. 

C ‘API* I j THE GRIGIONAL PRESSURE ( AT THE CENTER OF THE CLOUD). 

C ' A ZN ' IS THE HEIGHT OF THE CENTER OF THE CLOUD. 

C • A 0 M T • IS THE TOTAL PULLUTANT MASS IN EACH CLOUD. 

C 'PMAS‘» A VARIABLE DEPENDANT ON 'I' AS WELL AS 'L‘. IS THE MASS OF EACH 

C POLLUTANT IK EACH CLOUD. 

C L » 0 . C 

K * 0 

DO 1 1-1.20 
VOL ( I ) - 1 . 0 
PR OR ( 1 ) * 5 . 0 
TEMP(I )«1225.0 
ARMu ( I ) ■ 0. 0 
AIRMT(I)«U.O 
APMT ( I ) -0 . C 
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c j i j-i,n 
C J D< 1, Jl-Q.O 
Ct>£(l,J)»0.0 
Cl- Ft I > J ) * 0 . 0 
1 PINT-5 
CORTM-, FALSE. 

INTPS-.TRUE. 

GA8T.0-Q. 99 
QC-1. 1863E10 
0WATtP«3. 255789 
R£aD<5» JXM1) 

IPP-IPINT 
c K “6 . 2fc-5 
P-1.9&71 

Ul)«TF*wC-D*ATER 

RTS ( 1 ) *S QR T (9,8*JPHZ/T 1 1) ) 

ASPS(1)«.25*RTSU)*( <• 13822*3 <1)*RK*T(1>/P( 11/ Df>HZ)**.2 5) 
.1 /GA*TWO**Q» 75 
IF ( MAY* £0,1 ) GO TO 6 
DO 5 L»2,NLAf 
CALL TFIN0<TAHB»9ZKiUI 
CALL PFIND(PAM8j-8ZS(Li) 

PTS<L »-SQRT(9,8*OPHZ/TAH8) 

AZN(L)« (8ZN(L)+8ZN(L + U )/2.0 
CALL UFIND(UAMB» AZN(L) ) 

ASPS<L)*l»0/3»0*RTS(l)*((« 5529*QC*RK*TAHS/PAK8/DPHZ/UAMB > 
i *♦<1,0/3, on 

TFAR(L)«TFA+((dZN(L+l)**TFB)-(BZN(L)**TF3)) 

0<L >«OC*TFAR(L ) 

AVI < U»( <9ZN< L + 1>-3ZN( L H **3 , 0 ) *0. 19b3 5 

AZN(i)»(8ZN<2)“3ZN(l))/l,6A 

AVI < 1 ) * < < 4ZN< 1 )*GA1T«0> **3,0) *9. 1B8 79 

CNCMX«,01/AVK1) 

TFaRU)*TF 

T I f^E *— 0 T 

DO 7 L«1,NLAY 

CALL PFINO (PAM8, AZNCtn 

A P I < L ) -PAMS 

DO 7 1-1,6 
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PMASl*FRACA(I)*TFAR(L) 

IF (I.EQ.l.AND.L.EQ.l) PMASI-PMASI+5. 565E6 
APMT(L)»APMT(L)+PMASI 

7 PMAS ( L > I )«PMASI 

C *CPB E CPC' ARE CALCULATED SO THAT THE SPECIFIC HEAT OF POLLUTANT * I • 

C AT TEMPERATURE 'T* CORRESPONDING TO THE TEMPERATURE RANGE • J • IS EQUAL TO 
C *CPB(I#J) + CPC(I#J) * T*. 

C 'CPD# CPE# CPF* ARE CALCULATED TO FIND THE INTEGRAL OF 
C (MASS POLLUTANT + SPECIFIC HEAT) WITH RESPECT TO HEAT# OVER THE 
C TEMPERATURE RANGE FROM 298 K TO ANY TEMPERATURE UP TO 1300 K. 

C *CUMI ' IS THE SUMATIDN OF THE INTEGRALS OF THE PREVIOUS ( OR LOWER) 

C TEMPERATURE RANGES# ' T1 • IS THE BOTTOM TEMPERATURE OF THE RANGE# AND 
C *S1* IS THE SLOPE OF SPECIFIC HEAT/TEMPERATURE. 

C *CPA ( I# J ) * IS THE SPECIFIC HEAT OF PDLLUT ANT * I ' 

C IN TEMPERATURE RANGE *J*. 

DO 8 1*1# 6 

CUMI»-49.0*(2*CPA(I# 1)+.98*(CPA< I » 2 )-C PA ( I# 1 ) ) ) 

DO 8 J*l»ll 

T1 * 100. 0*J + 100.0 

C PA J*C PA ( I# J ) 

CPAJP*CPA( I# J+l) 

S1*(CPAJP-CPAJ)/100.0 

CPB(I#J)»CPAJ-S1*T1 

CPC(I#J)»S1 

IF U.NE.l) CUMI*CUMI+(CPAJ+CPA( I# J-l) )*50.0 
DO 8 L*l# NL AY 

CPD(L# J)»CPD(L#J)+(CUMI-T1*CPAJ+S1*T1+T1/2.0)*PMAS(L#I)/RMW(I) 

CPE CL# J)*CPE(L# J)+(CPAJ-T1*S1)*PMAS(L» I)/RMW(I ) 

8 CPF(L» J )*CPF(L# J)+(S1/2.0)*PMAS( L» I)/RMW(1) 

IF(INTPR)WRITE(6#212) UAMB# ( AS PS ( L ) #RTS ( L ) # L*l# NLA Y) 

212 FORMAT (* UAMB#ASPS#RTS ! *G15.6#/(1X#2G15.6)) 

A* ( ASPS (.1) *4. O/RTS (1)) **4.0 
TIME0-1-2.56E10/A 
IF (TIMED. LT.-1.0)TIME0*-1.0 
TIMEO-ACOS (TIMEO) 

ENTRY POINT 

C ABOUT TIME: THE TIME USED IN THE FOLLOWING SECTION OF THE PROGRAM IS 

C SET TO THE END OF THE INTERVAL# BUT THE TEMPERATURE RETURNED IS FOR THE 
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C BcGGINNING# AND THE CORRESPONDING TINE RETURNED AS • TL I H ' IN THE HAIN 
C PROGRAM MUST BE THE BcGGINNING. HENCE# • DT * IS SUBTRACTED FROM 'TIME' 

C BEFORE RETURN. 

C « N • IS THE TIME INTERVAL. 

C THIS PROGRAM USES A BINARY APROXIMATION ROUTINE# SETTING TD1 ( THE LOWER 
C RANGE) TO 200 K# AND T02 (THE UPPER) TO THE PREVIOUS TEMPERATURE (1275 
C ORIGIONALY ) . 1 AIRMQ' IS THE MASS OF AIR IN THE OLD CLOUD# 1 PMT • IS THE 

C POLLUTANT MASS# • VI » PI 1 ARE THE ORIGIONAL VOLLUME AND PRESSURE. «ZTE' 

C IS THE RISE IN THE CLOUD SINCE TIME 0. ‘DPR 6 PROR* ARE CALCULATED FROM 
C THE PRESSURE OF THE HEIGHT A ZN ( L ) # BEGINNING OF INTERVAL# CENTER OF CLOUD. 
15 TIME»TIMt+DT*2.0 
N*N+ 1 

DO 10 L * 1# NL AY 
TOWOO.O 
T02«TEMP(L)+50.0 
AIRMO-AR MO ( L ) 

PMT* APMT ( L ) 

VI-AVI(L) 

PI-API(L) 

ZTE-AZN(L)-BZN(L) 

CALL PFIND(PAMB»AZN(L) ) 

DPR(L)-PROR(L)-PAMB 
PROR ( L ) -PAMB 
TS«TIME*RTS(L) 

1 F ( C ORTM.AND.N.EQ. 1) TS=TIMEO 

C THE EQUATIONS FOR VOLUME AND VELOCITY ARE DIFFERENT FOR 'L»1‘. 

C 1 VOL £ D VOL • ARE CALCULATED FOR THE BEGINNING OF THE INTERVAL# 

C 'WN* IS THE VELOCITY# 'ZN« IS THE HEIGHT# AND 'VOBS' IS THE NEW VOLUME. 

IF (L.NE.l) GO TO 18 

VOBS- ( (ZTE + GAMTWO) *+3.0*4.18879 

DVOL (L)-VOBS-VOL(L) 

VOL(L)-VQBS 

WN«ASPS(1)*SIN(TS)*( (1-COS <TS ))*+(-. 75 )) 

ZN«AZN(1)+WN*DT 
AZN ( 1 ) -ZN 

V0BS-A.18879M(ZN*GAMTwD)**3.0) 

GO TO 19 

18 TT«BZN(L+1)-BZN(L) 

V0BS»TT*ZTE*ZTE*0. 78 54 
DVOL (L )«VQBS-VOL (L ) 
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VOL ( L )*VOBS 

WN“ASPS(L)+SlN(TS)*( ( 1-COS (TS))+*(-2. 0/3.0)) 

ZN*AZN(L)+WN*DT 
AZN(L)*ZN 
ZTE* ZN-B ZN ( L ) 

200 V0BS»TT+ZTE*ZTE*.785A 

19 CALL PFIND(PAMB*AZN<D) 

CALL TFIND(TAMB*AZN(L) ) 

C 'CPMIX' IS THE SPECIFIC HEAT OF THE CLOUD* 

C 'ZM' IS ONE OVER THE MOLECULAR WEIGHT* 

205 C ' A I RMP ' IS THE TOTAL HASS* 

C • Ml ' POINTS TO A TEMPERATURE RANGE EACH TIME THE FUNCTION IS CLACULATED. 

C 'Cl' IS THE MASS FRACTION OF A GIVEN POLLUTANT* 

C ' QL • <Q LOSS) IS CALCULATED IF 'N* NOT EQUAL TO 1. 

C » P VRK* CM1F* CZ* EC 1" ARE INTERMEDIATE VARIABLES INDEPENDANT OF 
210 C TEMPERATURE USED IN THE HEAT BALANCING EQUATION. 

C THE FIRST VALUES OF THE FUNCTION ARE CALCULATED* FOR 'TOl* AND 'T02'* 

C AND THEY ARE CHECKED TO MAKE SURE THEY ARE OF OPPOSITE SIGNS. IF NOT* 

C THE INTERVAL BETWEEN 'TOl* AND 'T02* IS DIVIDED INTO 50 SECTIONS* AND THE 
C VALUE OF THE FUNCTION IS PRINTED OUT AT EACH SECTION. IF THE VALUE EVER 
215 C CHANGES SIGN* THE PROGRAM BRANCHES ON* IF NOT* IT RETURNS. 

CPMIX«0.0 

ZM-0.0 

AIRMP-AIRMQ+PMT 
M1*INT(T02/100.0)-1 
220 DO 20 1*1*6 

CI*PMAS(L* Il/AIRNP 
ZM*ZM+CI/RMW< I ) 

IF(N.EQ.l) GO TO 20 

CPMIX*CPMIX+CI + (CPBU*M1) + CPC (I*H1)+T0 2) 

225 20 CONTINUE 

IF (N.EQ.l) GO TO 25 
C02*AIRM0+0. 21/AIR MP 
CN2-A1RM0+C.76/AIRMP 
ZM*ZM+C02/32.0+CN2/28.013 
230 CPMIX-CPMIX+AIRM0+.2A/AIRMP 

GAM-CPMIX/ (CPMIX-R+ZM) 

GAMR«1. O/GAM 

0L*IGAMR*VI*PI-GAMR*VI*(PI**GAMR)MPAMB+M1-GAMR) ) )*2.A21EA 
25 RKP»RK*ZM 
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PVRK-PAMB+VOBS/RKP 

CMlF-PVRKM-.24*TAMB + WN*WN/3 3 68.0) 
CZ*.24*(TAMB+AIRrtP+PVRK-AIRl1T{L))+QL-Q(L) 

C1--.244.PMT 
Ml ■ I NT (T01/100.0)-l 

TF1-CM1F/TD1 + CZ + CPD(L»M1> + T01*(C1 + CP£(L.M1) )+T014.T01*CPF(L.Ml) 
M1«INT(T02/100.0)-1 

TF2-CMlF/T02 + CZ + CPD(L»Ml)4-T02+{Cl + CPfc(L.Nl))+T024"T02*CPF(L.Ml) 
F1-AMAXHTF1.TF2) 

F2-AMINKTF1.TF2) 

1F(F1.GT. 0.0. AN0.F2.LT. 0.0) GO TO 70 
DO 90 NOIT-1.50 
TT"T01+NQIT*(TQ2-T01)/50.0 
M1«INT(TT/100.0)-1 

TF1-CM1F/TT+CZ+CPD(L#M1)+TT4.(C1 + CPE(L. Ml ) ) +TT4.TT+CPF (L. Ml) 
F1-AMAXHTF1.TF2) 

F2-AMINHTF1.TF2) 

WRITE(6.210)N0IT.F1.F2 
210 FORMATt* SSOJ FUNC :*I5, 2G15. 5 ) 

IF(F1.GT.0.0.AND.F2.LT.0.0) GO TO 92 
90 CONTINUE 
IF(INTPR) 

1 WRITE <6. 200) N0IT.CI11F.CZ»C1»CPD(L»M1 ).CPE(L. M1J.CPF (L. Ml). TOl. 

1 T02.TNEW.TF1.TF2.F1.F2.FNEW.RKP 
N*1 

TNEW-T02 

STOP 

92 T01-TT 

70 IF (TF1.GT.0.0) GO T071 

TNEW-T02 

T02-T01 

TQ1-TNEW 

C 'TOl L T02* ARE SET SO THAT 'TOl* CORRESPONDS TO THE POSITIVE VALUE AND 
C ' TQ2 * TO THE NEGATIVE (ABOVE). 

C IN THE MAIN ITERATION LOOP, THE FUNCTION IS CALCULATED FOR THE MIDDLE OF 
C THE INTERVAL. AND EITHER 'TOl* OR 'T02* IS SET TO THE MIDDLE OF THE 
C INTERVAL. DEPENDING ON THE FUNCTIONS* POLARITY. THE NEW MASS OF THE CLOUD 
C IS CALCULATED AND STORED. AND THE *TEMP» AND • DTK ' ARRAYS ARE FILLED. 

71 DO 75 NOIT-1.20 
TNEW-(T01*T02)/2.0 
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M1-INT(TNEW/100.0)-1 

FNEW«CM1F/TNEW+CZ+CPD(L,M1)+TNEW*!C1+CPE(L»N1) ) +TNEW*TNEW* 

1 C PF ( L / Ml ) 

IF (FNEW.GT.0.0) GO TO 72 
TQ2-TNEV 
GO TO 75 
72 T01-TNEW 

75 CONTINUE 

81 AIRN-PVRK/TNEW-AIRMP 

IF(N.GT.5.AN0.N.NE.IPR) GO TO 33 
IF (N.EQ.IPR.AND.L.EQ.NLAY) I PR -I PR + IPI NT 
IF ( INTPR ) 

1 WRITE (6j 200) NOIT#CMlFiCZ>Ci>CPD(L»Ml)>CPE(L*Ml)#CPF(L#Ml)>Taii 
1 T02,TNEW,TFl>TF2jFl,F2,FNEW,RKP 
200 FORMAT!* INTERMEDIATE OUTPUT: NOIT# CM IF, CZ> Cl > C P0» C PE ( L>M1) * 

1 /* CPF(LjM1)jT01j>T02#TNEW»TF1»TF2»F1jF2*/* FN£W*RKP*/ 

2 1X»I5»5G15.5»2(/1X»AG15.5)/1X»5G15.5) 

WRITE (6/ 100) N>L.TNEW>TIME>DT»AIRN>AIRMQ>NQIT»VOBSjZNjWN*QL*ZMj 
1 CPMIX>GAM,0<L)»PAMB# TAMB > PMT# P I> VI, AIRMT(L) 

33 AIRMT!L>«AIRMTIL)+AIRN*TAMB 
ARM01L)-AIRMQ+AIRN 
DTK!L)»TEMP!L J-TNEW 
TEMP ( L ) -TNE W 
10 CONTINUE 

T1ME-TIME-DT 
IF(N.EQ.l) GO TO 15 

100 F0RMAT!1X,12(AH****)//* TIME INTERVAL * *15/* LAYER -*15/ 

1 * NEW TEMPERATURE « *G15.5/* TIME « +G15.5/ 

2 * DELTA TIME ■ *G15.5/* MASS OF NEW AIR INGESTED INTO CLOUD - * 

3 G15.5/* MASS OF OLD AIR IN CLOUD - *G15.5/ 

A * NO. OF ITERATIONS « *15/* OBSERVED VOLUME - *G15.5/ 

5 * HEIGHT OF CLOUD » *G15.5/* CLOUD VELOCITY - *G15.5/ 

6 * OL - *G15.6/+ ADDITIONAL VARIABLES* ZM, CPMI X , GAM, Q 1 L ) , P AMB, 

7 TAMB, PMT»PI,VI,AIRMT (L)*2(/1X» 5G15.6 ) ) 

110 FORMAT! / / 1 X, 70 ( 1H* ) / /* NO CONVERGENCE */ /IX# 70 ( 1H* ) ) 

60 RETURN 

END 



SUBROUTINE TFIND ( TF> ZA ) 

COMKON/CB/ATZ, ATF 
DIMENSION ATZ(100),ATF<100) 

DO 10 1-1,100 
IP-I+1 
ZHI-ATZUP > 

IF(ZA.LT.ZHI) GO TO 12 

CONTINUE 

TLOW-ATF (I ) 

THI-ATFUP ) 

ZLOW-ATZ (I ) 

TF«TLOW+(ZA-ZLOW)*(THI-TLOW)/(ZHI-ZLOW ) 

RETURN 

END 


SUBROUTINE PFIND(TF,ZA) 
COMMON/CC/APZ,APF 
DIMENSION APZ(100),APF( 100 ) 

DO 10 1-1,100 

IP-I+1 

ZHI-APZ(IP) 

IF(ZA.LT.ZHI) GO TO 12 

CONTINUE 

TLOW-APF ( I ) 

THI-APF (IP) 

ZLOW-APZ ( I ) 

TF-TLCW+IZA-ZLOW )*(THI-TLOW) / (ZHI-ZLOW ) 

RETURN 

END 

SUBROUTINE UF IND ( TF, Z A ) 

COMMQN/CD/AUZ, AUF 
DIMENSION AUZ( 100), AUF ( 100) 

DO 10 1-1,100 
IP-I+1 
ZHI-AUZIIP ) 

IF(ZA.LT.ZHI) GO TO 12 

CONTINUE 

TLOW-AUF ( I ) 

THI-AUF(IP) 

ZLOW-AUZ ( I ) 

TF-TLQW+<ZA-ZLOw )+<THI-TLOW> / (ZHI-ZLOW ) 

RETURN 

END 
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SHAPE 


CENTER LOCATION 


VOL 

X 

y 

z r 

(meters) 

AD 

1 

33.12 

-46.4 

96.5 

533.9 

193 

2 

42.88 

-73.1 

209.8 

533.9 

33.6 

3 

67.45 

-153.7 

265.7 

533.9 

78.4 

4 

174.38 

-553.1 

431 

533.9 

252.1 

5 

197.9 

-685.6 

586.9 

533.9 

59.8 

6 

208.8 

-807.5 

639.1 

533.9 

44.6 




UPPER AIR METEOROLOGICAL DATA 
Dec. 10, 1974 0631Z 


BOUNDARY LAYER 
meters 

WIND DIRECTION 
degrees 

WIND SPEED 
met/s 

TEMPERATURE 
degrees C 

PRESSURE 

mb 

RELATIVE HUMIDITY 
percent 

4.88 

310.0 

3.08 

7.9 

1023.0 

62.0 

192.99 

339.0 

9.77 

9.1 

1000.0 

66.0 

226.52 

341.0 

9.77 

9.5 

996.0 

67.0 

304.88 

345.0 

10.80 

8.9 

986.6 

68.0 

557.01 

345.0 

10.28 

6.7 

957.0 

60.0 

616.77 

355.0 

10.28 

7.8 

950.0 

57.0 

661.28 

354.0 

10.28 

8.6 

945.0 

55.0 

914.63 

346.0 

10.28 

7.9 

916.5 

40.0 

1063.11 

340.0 

11.31 

7.4 

900.0 

40.0 

1219.51 

332.0 

11.83 

7.1 

883.2 

39.0 

1524.39 

314.0 

12.85 

6.5 

851.1 

37.0 

1533.23 

314.0 

12.85 

6.4 

850.0 

37.0 

1829.27 

306.0 

13.88 

5.8 

820.1 

37.0 

1870.73 

305.0 

13.88 

5.8 

816.0 

37.0 

2028.96 

318.0 

13.88 

4.3 

800.0 

37.0 

2134.15 

319.0 

13.88 

3.6 

790.1 

37.0 

2238.72 

310.0 

15.42 

3.0 

780.0 

37.0 


Stability Category: Stable Air 

Brunt Voisala Frequency = 0.0172 (sec -1 ) 

Mixing Layer Depth: 665 (meters) 
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SOURCE DATA 
August 20, 1975 
2122 Z KSC 


SHAPE 

CENTER LOCATION r AD POLLUTANT STRENGTH IN EACH VOLUME 


VOL 

X 

y 

z 

(meters) 


CO 

C0 2 HC1 

(milligrams) 

A1 2 0s 

1 

-25.6 

5.7 

83.3 

831.4 

166.5 

4.0331436E+07 

4.1927008E+06 

3.0401153E+07 

4.3937757E+07 

2 

-43.2 

6.6 

193.3 

831.4 

53.6 

2.2114947E+07 

2.2989848E+06 

1.6669872E+07 

2.4092402E+07 

3 

-78.5 

6.6 

262.5 

831.4 

84.8 

5. 7781849E+07 

6.0067785E+06 

4.3554979E+07 

6.2948537E+07 

4 

-237.2 

-4.6 

457.4 

831.4 

304.9 

7.6493345E+08 

7 . 9519535E+07 

5.7659387E+08 

8.3333160E+08 

5 

-252.7 

-6.6 

624.1 

831.4 

28.6 

1.5260289E+08 

1.5864009E+07 

1.1502948E+08 

1.6624821E+08 

6 

-377.2 

-18.5 

776.5 

831.4 

276.2 

2. 7565466E+09 

2 . 8655998E+08 

2.0778381E+09 

3 . 0030291E+09 

7 

-425.1 

-23.8 

1002.6 

831.4 

176 

3.1968884E+09 

3.3233620E+08 

2.4097604E+09 

3.4827450E+09 

8 

-461.3 

-21 

1155.1 

831.4 

128.9 

2. 9295459E+09 

3 . 0454432E+08 

2.2082421E+09 

3. 191497 5E+09 

9 

-529.9 

13.9 

1374 

784.8 

308.9 

9 . 0534840E+09 

9.4116533E+08 

6.8243627E+09 

9.8630206E+09 

10 

-131.7 

594 

1553.8 

668.7 

58.8 

1. 6301156E+09 

1 . 6946054E+08 

1.2287535E+09 

1.7758758E+09 

11 

-167.7 

584.4 

1706.3 

570.8 

246.1 

5.3512145E+09 

5.5629164E+08 

4.0336548E+09 

5.8297048E+09 

12 

-243.9 

879.6 

1908.6 

441.4 

158.5 

2.2781997E+09 

2.3683287E+08 

1. 7172683E+09 

2.4819098E+09 
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UPPER AIR METEOROLOGICAL DATA 
Aug. 20, 1975 2046Z 


BOUNDARY LAYER 
meters 

WIND DIRECTION 
degrees 

WIND SPEED 
met/s 

TEMPERATURE 
degrees C 

PRESSURE 

mb 

RELATIVE HUMIDITY 
percent 

4.88 

110.0 

3.60 

28.7 

1018.3 

75.0 

166.50 

95.0 

4.11 

26.8 

1000.0 

84.0 

220.12 

91.0 

4.11 

26.2 

994.0 

83.0 

304.90 

89.0 

4.11 

25.5 

984.5 

83.0 

609.80 

83.0 

2.60 

24.5 

951.1 

62.0 

638.41 

82.0 

2.60 

24.5 

948.0 

62.0 

914.60 

87.0 

1.03 

21.8 

918.5 

69.0 

1090.60 

80.0 

.51 

20.6 

900.0 

74.0 

1219.50 

109.0 

.51 

20.0 

886.9 

69.0 

1524.40 

169.0 

1.54 

17.7 

856.1 

68.0 

1583.20 

166.0 

1.54 

17.4 

850.0 

69.0 

1829.30 

162.0 

1.54 

16.3 

826.1 

70.0 

1987.80 

167.0 

3.08 

15.4 

811.0 

70.0 

2100.30 

170.0 

4.11 

15.5 

800.0 

71.0 

2211.30 

171.0 

5.66 

15.6 

790.0 

73.0 

2439.00 

159.0 

4.63 

11.7 

769.0 

77.0 


Stability Category: Stable Air 

Brunt Voisala Frequency = 0.008 (sec -1 ) 

Mixing Layer Depth: 1988 (meters) 


SOURCE DATA 
March 14, 1976 
0127 Z March 15 KSC 

SHAPE 

CENTER LOCATION POLLUTANT STRENGTH IN EACH VOLUME 


VOL 

X 

y 

z 

(meters) 

r 

AD 

CO 

CO a HQ1 

(milligrams) 

AI2O3 

1 

-10.5 

-29.0 

91.5 

608.7 

182.3 

6.74629 E+07 

7.0123 E+06 

5.07786 E+07 

7.3508 E+07 

2 

-64.3 

-73.9 

244.2 

608.7 

122.6 

1.61505 E+08 

1.6787 E+07 

1.215631E+08 

1.75977E+08 

3 

-336.9 

-154.9 

457.4 

608.7 

304.9 

2.095226E+09 

2.17783E+08 

1.577052E+09 

2.28297E+09 

4 

-354.6 

-155.7 

617.4 

608.7 

15.2 

2.186213E+08 

2.27241E+07 

1.645537E+08 

2.38211E+08 

5 

-831.6 

-105.1 

769.8 

608.7 

289.6 

6. 607105E+09 

6.86790E+08 

4.973090E+09 

7.19914E+09 

6 

-987.1 

-56.9 

1001.4 

576. 

173.5 

6. 203409E+09 

6.44799E+08 

4.669233E+09 

6. 75927E+09 



UPPER AIR METEOROLOGICAL DATA 
March 14, 1976 

0127Z, March 15, 1976 (2027 EDT, March 14, 1976) 


BOUNDARY LAYER 
meters 

WIND DIRECTION 
degrees 

WIND SPEED 
met/s 

TEMPERATURE 
degrees C 

PRESSURE 

mb 

RELATIVE HUMIDITY 
percent 

4.88 

0.0 

1.54 

19.9 

1020.7 

87.0 

182.30 

40.0 

6.17 

19.9 

1000.0 

98.0 

304.88 

60.0 

5.66 

19.6 

985.9 

97.0 

609.76 

87.0 

5.66 

18.2 

951.8 

99.0 

625.00 

87.0 

5.15 

18.1 

950.0 

99.0 

914.63 

105.0 

3.60 

16.9 

918.6 

84.0 

1088.10 

130.0 

3,08 

15.3 

900.0 

98.0 

1219.50 

154.0 

3.08 

14.9 

886.4 

97.0 

1524.40 

211.0 

3.60 

13.1 

855.1 

100.0 

1572.60 

218.0 

4.12 

12.8 

850.0 

100.0 

1829.30 

233.0 

5.66 

11.5 

824.7 

95.0 

2081.10 

239.0 

6.70 

10.0 

800.0 

94.0 

2134.20 

240.0 

6.70 

9.6 

795.2 

95.0 

2439.00 

248.0 

7.20 

7.9 

766.6 

77.0 

2615.50 

252.0 

7.20 

6.2 

750.0 

79.0 


Stability Category: Stable Atmosphere 

Brunt Voisala Frequency = S 

_ /g_ 3<K 1/2 
Xo 3z ; 

0.014 (sec -1 ) 


Mixing Layer Depth: 1088 (meters) 



z 



/ 

/ 

FIGURE C-l DESCRIPTION OF THE SOURCE DATA 

The stabilized cloud Is layered Into several 
subclouds, each of which resembles a pancake. 
The center (x,y,z), the height AD and the 
radius r of each subcloud are given by means 
of referring to a cartesian coordinated sys- 
tem set at the launch pad. The total amount 
of each pollutant species in each subcloud is 
given in the source data tables. It has been 
assumed that the concentration distribution 
in each subcloud is Gaussian in the horizontal 
plane, xy, and uniform in the vertical z 
direction. 
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APPENDIX D 


D,i 


C. 

i 


v 

_ i 


D. 

x 

d_ 

m 


ij 


w 


K 


gas 


K. . 
ij 

K?. 

ij 

K , K , 

x y 

K„ 


B 


L. . 
1J 



SYMBOLS 

mean concentration of contaminant 
drag coefficient for particle size class i 
mass fraction of species i in cloud 
specific heat at constant pressure 
specific heat at constant density 
fluctuating concentration 

particle diffusion coefficient (cm 2 sec -1 ) 
mass average diameter (cm) 
saturation vapor pressure (mb) 

Coriolis parameter (sec -1 ) 

mean thermal speed of ith and jth sized particles (cm sec -1 ) 
gravity constant (cm sec -2 ) 
thickness of flow field 
cloud thickness 

effective heat release per unit mass of liquid propellant (cal g -1 ) 
effective heat release per unit mass of solid propellant (cal g -1 ) 
latent heat of water vaporization at room temperature (cal g -1 ) 
gas constant (atm»m 3 «K -1 ) 

eddy diffusivity tensor (m 2 sec -1 ) 
collision kernel 

eddy dif fusivities (m 2 sec -1 ) 

Von Karman constant 
Boltzmann's constant (ergs K"* 1 ) 

Monin-Obukhov length scale 

rate of collision between particles of component A/B, size i, with 
particles of A/B, size j 

rate of collision of particles of component B, size i, with particles 
of component B, size j 
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4 

M 

s 

M 

w 

MW. 


M 


total 


m 


m . 
air 

in o 

N' (d) 
n 
n 0 
P 

^amb 

Q 

4 



R 

g 

RH 


R. . 
ij 

S 


s 

T 

T22 > T3 3 



t 


fc f 

4 

u 

p 
u I 


rate of collision of component A/B, size i, with particles of 
component B, size j 
mixing length (m) 

mass flow rate of liquid engine (g sec -1 ) 
mass flow rate of solid engine (g sec -1 ) 
mass flow rate of poured water (g sec -1 ) 
molecular weight of species i in cloud 
actual initial mass of aerosol (g) (see eq. 14) 
mass (g.) 

mass rate of air entrained due to turbulent mixing (g sec -1 ) 

initial average particle mass (g) 

initial particle number distribution 

number of moles in a unit volume 

initial (total) particle number density 

pressure in cloud (atm) 

ambient environmental pressure (atm) 

source strength (ppm) 

effective available heat (total effective heat content in cloud) 
(cal) 

particle Reynolds number 
flux Richardson number 
gradient Richardson number 
relative humidity (%) 

Lagrangian correlation tensor 
constant mean shear gradient (sec -1 ) 

Brunt Vaisala frequency (sec -1 ) 
absolute temperature (K) 

Lagrangian integral time scales in crosswind and vertical direction 
respectively (sec) 

absolute temperature of ambient air near ground (K) 
time 

firing time during which cloud is formed (sec) 
time for initial cloud rise speed (sec) (see eq. 5) 
pseudovelocity (m sec -1 ) 

fluctuating velocity component in i direction (m sec 1 ) 
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V 

S,1 

W 

W m 

w'q' 

X 0 

x 0 

Xoo 

z 0 

Z. 

l 

Z 

m 

Z(t) 

z 

z 0 

a 


y 

Y' 


Y 

Y 

Y 

Y 

Y 
6 


d 

g 

P 

T 

u 


e 

0oo 

X 

y 


V 

p 

p 


friction velocity (m sec -1 ) 
wind speed (m sec -1 ) 
geostrophic wind speed (m sec -1 ) 
observed cloud volume (m 3 ) 

mixed Eulerian/Lagrangian fluctuation velocity (m sec -1 ) 
sedimentation velocity of i size particle (cm sec -1 ) 
cloud rise speed (m sec -1 ) 
mixing ratio 
heat flux 

alongwind distance from pollutant point source 

point source location 

fictitious point source location 

initial cloud center location 

thickness of the PBL 

predicted stabilization height 

cloud altitude at time t 

height of source location 

roughness length 

mass class size 

entrainment constant 

specific heat ratio 

adiabatic lapse rate (K m -1 ) 

constant for absorption strength of lower boundary 

absorption factor 

eddy conductivity (m sec -1 ) 

constant for absorption strength of upper boundary 

mean free path for particles (cm) 

rate of energy dissipation (m 1 * sec -1 ) 

total mass of pollutants at altitudes z and time t 

molecular mean free path (cm) 

viscosity of air (g cm sec -1 ) 

eddy viscosity (m 2 sec -1 ) 

set of post-afterburning chemical species 

particle density (g m -3 ) 
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axr 


g 

tfo 


ij 


w 


o , o 
x y 


x 0 


Yo 


z 0 


density of ambient air near ground (g m 3 ) 
density of air (g ml -1 ) 
initial standard deviation 
standard deviation of wind azimuth angle 

standard deviation of wind azimuth angle at a reference height near 
ground 

standard deviation of wind elevation angle 

standard deviation of wind elevation angle at a reference height near 
ground 

dispersion tensor 

standard deviations of turbulent wind fluctuation components in cross- 
wind direction 

standard deviations of turbulent wind fluctuation components in 
vertical direction 

standard deviations of mean concentration dispersions in x, y, z 
directions 

standard deviation of initial pollutant distribution in alongwind 
direction 

standard deviation of initial pollutant distribution in crosswind 
direction 

initial standard deviation of initial pollutant distribution in 
vertical direction 
time step size (sec) 
nondimens ional wind shear 

vertical gradient of virtual potential temperature 


T 

<f> 

/ Az , 

9<f>/3z 

Subscripts : 
air ambient air 

amb ambient 

obs observed 

max maximum 
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TABLE I 


POST-AFTERBURNING PLUME COMPOSITION 


AlaOa 

n 2 

HC1 

H 2 0 

C0 2 

Ffi 2 O 3 

Mass Frac . 0 . 1214 

0.4882 

0.0867 

0.1344 

0.1677 

0.0016 


Y. 

l 


TABLE II 



COMPARISON OF 

AFTERBURNING 1 

COMPOSITION 


OF 

MAJOR SPECIES 



Complete 

Titan III a 

Shuttle 


burning 

0.9 km * 5 

0.7 kmb 


( 8 ) 

(g) 

( 8 ) 

HC1 

21.6 

18.7 

19.0 

CO 2 

41.8 

41.4 

41.4 

H 2 0 

33.5 

29.0 

28.7 


Value at 1 km downstream of the nozzle exit plane. 

Altitude of the rocket motor; determines the ambient undisturbed 
atmospheric condition and the speed of the rocket. 
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TABLE III 


THE RANGE OF MULTILAYER DIFFUSION MODEL PARAMETER VARIATIONS 
IN STEWART AND GROSE STUDY 


Lateral diffusion exponent, a 

Vertical diffusion exponent, 3 

Source strength in each layer, 

Q k , ppm-m 2 

Source-strength distribution in 

each layer, Q , ppm-m 2 
K 

Mixing-layer division (layer grid 
spacing) , K 

Stabilized cloud geometry 

Specific energy release from rocket 
motor, h g , cal/gm 

Meteorology 

Mixing-layer depth, , meters 


0.5 to 1.5 
0.5 to 1.5 
10 3 to 2 x 10 7 

Uniform; Gaussian; step function 


7 layers to 28 layers 


Right circular cones; right circular 
cylinders: both with equal volumes, 
with and without equal mass loadings 

690 to 2980 


Low-level sea-breeze regime; fall 
fair-weather regime at Kennedy Space 
Center 

300 to 750 
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TABLE IV 


STANDARD DEVIATIONS OF WIND AZIMUTH ANGLE cr A 

A 0 

AND ELEVATION ANGLE a_ AT REFERENCE HEIGHT 

Eo 

4 m ABOVE GROUND 

10 December 1974 20 August 1975 

(Towers 110, 313) (Towers 110, 308, 311, 313) 

Source CT A 0 g E 0 g A 0 °E 0 

i.l 2 2 10 10 

i. 2 8 4 11 7 

i. 3 8 5 11 8 

1.4 88 23 23 


TABLE V 

PEAK OF MAXIMUM CONCENTRATION (ppm) AT GROUND FROM MDM 
BY VARYING THE STANDARD DEVIATION OF THE INITIAL POLLUTANT 


DISTRIBUTION IN DECEMBER LAUNCH CASE 3 
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TABLE VI 

BRIEF DESCRIPTION OF SELECTED 
DIFFUSION MODELS 

Models (Sponsor /User) Modeling Technique Turbulence Schemes 

a 

ADPIC (ERDA/LLL) Semi-Trajectory Particle Similarity Method for 

Diffusion Model Diagonal Eddy Diffusivity 

DISF (NASA/AEROCHEM) Analytical Model for Lagrangian Method 

Uniform Shear Flow Eddy Diffusivity Tensor 

MDM (NASA/AEROCHEM) Gaussian Plume Model Cramer’s Constants^ 

Ignoring X-Dif fusion 

METS (VAFB/SAI) Gaussian Plume Model Turner -Pas quill^ 

Ignoring X-Diffusion 

TREATS (ERDA/NUS) Gaussian Distribution in Moment Method for a x , Oy 

Each Horizontal Plane P-G or Blackadar Type for o z 

Can be changed for any specified diagonal eddy diffusivity. 

^Alongwind dispersion a x is based on Tyldesley and Wallington's work. 
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TABLE VII 


RATIO OF TURBULENCE PARAMETERS 


0 

u 

a 

V 

a 

w 

u* 

Reference 

Description 

1.68 

1.28 

1.19 

1 

Champagne 
et al- 

Homogeneous shear flows in 
a wind tunnel 

2.0 

1.7 

1.25 

1 

Cramer 5 9 

Atmospheric surface layer 

1.89 

1.32 

0.97 

1 

Klebanoff 60 

Wind tunnel boundary layer 
on a smooth surface 

2.0 

1.4 

1.27 

1 

Wyngaard 
. .6 1' 
et al 

Atmospheric surface layer, 
modified data 

1.82 

1.3 

1.12 

1 

So and 
Mellor 62 

Surface layer in wind tunnel 
smooth wall 

2.4 

1.43 

1.03 

1 

Comte-Bellot 6 3 

Turbulent flow between 
parallel plates 

2.2 

1.7 

0.9 

1 

. 6 4 

Hinze 

Turbulent pipe flow 

2.5 

1.7 

1.25 

1 

4 3 

Yaglom 

Atmospheric surface layer 

2.5 

2.2 


1 

Blackadar 
et al 31 

Atmospheric boundary layer 
at KSC 
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TABLE VIII 


COMPARISON OF HC1 CONCENTRATION AMONG MODELS AND OBSERVATIONS AT 
MEASUREMENT LOCATION P-10 FOR TITAN III LAUNCH OF AUGUST 20, 1975 
P-10 

Measurement ADPIC a DISF METS b MDM TREATS 3, 


Maximum 

instantaneous 

concentration 

ppm 


1.4 x 10" 2 0,02 x 10 -2 1.7 x 10- 2 


54. x 10" 2 0.68 x 10" 2 


Integrated 

concentration 7 >8 39 >18 241 1 

(dosage) 
ppm- sec 


a Value obtained after justification of coordinate system for comparison (see text). 
Walue guessed from available data. 


TABLE IX 

COMPARISONS OF THE DILUTION RATIOS OF THE UPPER LEVEL CONCENTRATIONS 
FOR MODEL PREDICTIONS AND AIRBORNE MEASUREMENT FOR 10 DECEMBER 1974 


Time 

ADPIC 

DISF 

MDM 

TREATS 

Measurement 

Earlier 

0.81 

0.94 

< 0.11 

0,60 

— 

Later 

0,83 

0.89 

0,51 

0.49 

0,60 to 0.82 
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TABLE X 

COMPARISON OF HC1 CONCENTRATION AMONG MODELS AND OBSERVATIONS AT MONITORING SITES 

P-2, P-3, P-4 FOR DECEMBER CASE 


MDM 




Measurement. 

ADPIC a 

DISF 

b 

c 

d 

METS e 

TREATS 


Maximum 









P-2 

instantaneous 

concentration 

0.35 

0.03 

0.07 

0.12 

1.14 

0.50 

T-* ~ 

0.11 


Integrated 

concentration 

19.5 

> 6 

15 

27 

240 

62 

25 

58 


Maximum 









P-3 

instantaneous 

concentration 

0.022 

0.01 

0.013 

0.14 

0.55 

0.52 


0.03 


Integrated 

concentration 

6,2 

S 6 

5 

45 

183 

83 

< 0.7 

7.5 


Maximum 









P-4 

instantaneous 

concentration 

0.50 

0.075 

0.11 

0.28 

1.26 

0.33 


0.21 


Integrated 

concentration 

15.2 

> 45 

15 

33 

145 

28 

< 40 

56 


a The integrated concentration was based on the first 55 min integration. 
^Vertically varied o ^ . 



^alue was guessed from limited results. 



TABLE XI 

COMPARISON OF HC1 CONCENTRATION USING OBSERVED CLOUD 
AS INPUT AT MONITORING SITE P-10 FOR TITAN III LAUNCH 

OF AUGUST 20 , 1975 

MDM 

Measurement DISF a DISF^ Model 3 Model 4 

Maximum 

instantaneous 1.4 x 10 -2 1.5 x 10 -2 0.27 x 10~ 2 63. x 10 -2 63. x 10~ 2 

concentration 
ppm 

Integrated 

concentration 7 33 5.5 117 114 

(dosage) 
ppm-sec 


3 

After adjustment, 
b 

Before adjustment. 
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TABLE XII 

COMPARISON OF NUMERICAL AND ANALYTICAL SOLUTIONS 
TO THE COAGULATION RATE EQUATIONS 



time 

N 

total 

d N,90 

^,50 

d N,10 

d_ 

m 


(sec) 

(ml' 

- 1 ) 

<ym) 

(Uni) 

(ym) 

(ym) 

Analytical 


0 

1.00 x 

10 6 

6.44 

9.43 

13.9 

10.0 

Solution 

1.25 

x 10 5 

2.60 x 

10“ 

16.4 

29.9 

44.6 

33.8 


1.7 

x 10° 

1.96 x 

10 3 

37.9 

70.8 

105.5 

80.0 

Numerical 









Solutions 









a = 2, T = 0.01 

1.25 

x 10= 

2.59 x 

10“ 

15.5 

29.0 

44.5 

33.8 

a = 2 , x = 0. 1 

1.25 

x 10 s 

2.56 x 

10“ 

16.0 

29.6 

44.4 

33.6 


1.7 

x 10 6 

2.00 x 

10 3 

35.5 

67.0 

103.0 

79.5 

a = 10, x = 0.1 

1.25 

x 10 5 

2.56 x 

10“ 

13.0 

25.1 

39.0 

33.6 


1.7 

x 10 6 

2.04 x 

10 3 

26.5 

45.8 

89.0 

78.7 


TABLE X2II 

VALUES FOR THE COAGULATION KERNEL, K^ 
(10 -1 ° ml particle -1 sec -1 ) 


Particle 3 
diam (ym) 

0.010 

0.026 

0.089 

0.30 

1.0 

3.4 

0.10 

19.5 

43.4 

219.2 

940.0 

3351.0 

11390.0 

0.026 

43.4 

24.3 

49.3 

155.5 

515,0 

1725.0 

0.089 

219.2 

49.3 

16.1 

22.6 

60.0 

189.9 

0.30 

940.0 

155.5 

22.6 

9.06 

12.7 

39.2 

1.0 

3351.0 

515.0 

60.0 

12.7 

6.91 

35.4 

3.4 

11390.0 

1725.0 

189.9 

39.2 

35.4 

6.28 


a Particles have unit density. 
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TABLE XIV 

INITIAL ALUMINA AND DEBRIS PARTICLE DISTRIBUTIONS 


Total , , 

Particles N,50 m 


Case 

(ml - 

- 1 ) 


Distribution, particles ml -1 ym -1 

(ym) 

(ym) 

Alumina 

Dawborn 

5.6 x 

10 7 

2.7 x 10 9 [d 2 exp (-10. 5 d 1 / 2 ) + 1 x 

10“ d 3 exp(-A1.25d)] , d ^ 0.02 ym 

0.088 

0.16 

Kreaulte 

A. 8 x 

10= 

6. A x 10 s [d 3 exp (-3.0 d) + 1.0A x 

10 s d“ exp(-25.39d l/a )], d^0.02 ym 

0.16 

0.77 

Varsi 

8.8 x 

10 7 

1.10 x 10 7 x < 

f d“ l - 2S 

3.6 x 10 -2 d -3 * 2 

0.02 ym ^ d ^ 0.06 ym 
0.06 ym ^ d ^ 0.10 ym 

0.066 

0.12 

Debris 

With Dawborn 

8.6 x 

10 7 


^ 1.58 x 10- 3 d-“- 8 
1.05 x 10 6 ' d -3 * 27 

0.10 ym ^ d ^ 10 ym 

0.1 ym ^ d ^ 1000 ym 

0.1A 

1.36 

(1000 x alumina) 
With Kreautle 

3.6 x 

10 7 


3.0 x 10 s d -3 -'* 7 

0.1 ym ^ d ^ 1000 ym 

0.13 

0.8A 

(100 x alumina) 
With Kreautle 

A. 9 x 

10 7 


6.8 x 10 s d” 3 - 20 

0.1 ym ^ d ^ 1000 ym 

0.1A 

1,6 A 

(1000 x alumina) 
Varsi 

1.6 x 

10 8 


1.6 ^ 10 s d" 3 - 33 

0.1 ym ^ d ^ 1000 ym 

0.13 

1.11 


TABLE XV 

ALUMINA LOADINGS BELOW THE CENTER 
OF THE GROUND CLOUD 


Distribution 
Dawborn, debris 
1000 x alumina 

Kreautle, debris 
100 x alumina 

Kreautle, debris 
1000 x alumina 

Varsi, debris 
1000 x alumina 



0.9 


4.5 

1.4 
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TIME, min 


FIGURE 1 PHYSICAL 


ENTRAINED AIR MASS RATE 

gm tec-* TEMPERATURE 
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DURING CLOUD RISE 
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VIRTUAL POTENTIAL 
TEMPERATURE, K 



rad. 



WIND SPEED, 

m mc " 1 



RELATIVE HUMIDITY, 
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FIGURE 2 VERTICAL PROFILES OF TEMPERATURE, WIND SPEED, WIND DIRECTION, AND RELATIVE 
HUMIDITY AT KENNEDY SPACE CENTER FROM RAWINSONDE SONDING MEASUREMENT 

FOR 10 DECEMBER 1974, 02:31 EDT 




HEIGHT 





MAXIMUM CONCENTRATION, ppm 
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FIGURE 4 MAXIMUM GROUND LEVEL CONCENTRATION PREDICTED BY MDM MODEL 4 
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RATIO OF MAXIMUM CONCENTRATION, C/C(y p «0) 
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GROUND LEVEL HCI MAXIMUM DOSAGE , ppm-ttc 



FIGURE 6 VARIATIONS OF PREDICTED GROUND LEVEL HCl MAXIMUM DOSAGE DUE TO VARIATIONS OF 

SURFACE ABSORPTION FACTOR y AND DIFFUSION PARAMETERS a AND cf FOR THE 

p A E 

10 DECEMBER 1974 LAUNCH 


Ln 






HEIGHT 






RATIO OF SOURCE STRENGTH, Q(y)/Q(Oj 64 ) 




ENTRAINMENT CONSTANT RATIO, y/y Q {y Q - 0.64) 


FIGURE 8 VARIATIONS OF PREDICTED SOURCE STRENGTH IN 
CLOUD USING DIFFERENT ENTRAINMENT CONSTANTS y FOR 
THE 10 DECEMBER 1974 AND 20 AUGUST 1975 CASES 
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FIGURE 9 MAP OF KSC AREA 


L5 


Tower heights are 16.5 m except tower 
110 (62 m) and tower 313 (150 m) 




z/Zi 



FIGURE 10 MODELED REYNOLDS STRESS IN THE PBL FOR 
THE 10 DECEMBER 1974 AND 20 AUGUST 1975 CASES 


Dashed line is the actual input between 0.2 and 
0.5 z/Z. for the DISF calculation 
for the August case 
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FIGURE 11 VERTICAL PROFILES OF TEMPERATURE, WIND SPEED, WIND DIRECTION, AND RELATIVE 
HUMIDITY AT KENNEDY SPACE CENTER FROM RAWINSONDE SONDING MEASUREMENT 

FOR 14 MARCH 1976, 20:27 EDT 














GROUND LEVEL HCI MAXIMUM DOSAGE, ppm-»tc 



DISTANCE, km 


FIGURE 14 COMPARISONS OF GROUND LEVEL HCI MAXIMUM DOSAGE 
PREDICTIONS FOR THE 20 AUGUST 1975 CASE 

MDM.l refers to results using measured a ^ ; 

O 

MDM.2 to from theoretical interpolation; 

MDM.3 to a given in Ref. 23 

n. 
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MAXIMUM CONCENTRATION, ppm 



FIGURE 16 COMPARISONS OF MAXIMUM GROUND LEVEL CONCENTRATION 
PREDICTIONS FOR THE 10 DECEMBER 1974 CASE 

MDM Model 3.1 and MDM Model 4.1 calculated using input 
from SAI preprocessor version. MDM ii calculated 
using monotonically decaying a and a E « 
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DISTANCE) km 

FIGURE 17 COMPARISONS OF GROUND LEVEL HC1 MAXIMUM DOSAGE 
PREDICTIONS FOR THE 10 DECEMBER 1974 CASE 

MDM 1.3 calculated using a. = 8, a = 5. MDM ii calculated using 

A Hi 

O O 

monotonically decaying a^and cr^. DISF 1 and DISF 2 calculated 
including and ignoring variations in wind direction, respectively. 
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CONC. 



3 5 7 9 II 13 

DISTANCE, Km 


FIGURE 18 COMPARISONS OF MAXIMUM GROUND LEVEL CONCENTRATION 
PREDICTIONS FOR THE 14 MARCH 1976 CASE 

MDM 1 calculated using interpolated a ^ . 

O 

MDM 2 calculated using measured . 

O 


O' 

LO 









Each line represents one minute. Numbers in parentheses 
represent aircraft pass numbers. 
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FIGURE 21 AIRCRAFT ALTITUDE ON 20 AUGUST 1975 FOR CLOUD PASSES 1 THROUGH 15 

Each pass of one minute duration. Numbers represent 
aircraft pass numbers. 
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TIME, sec. 

FIGURE 22 COMPARISONS OF AIRBORNE SAMPLING DATA WITH (a) THE TREATS RESULTS AND 
(b) THE ADPIC RESULTS FOR THE 20 AUGUST 1975 CASE 

Key: Model Upwind relocation 

Measured data Downwind relocation 








CONCENTRATION, ppm 


Warn 


10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60 

TIME, sec. 


FIGURE 23 COMPARISONS OF AIRBORNE SAMPLING DATA WITH RESULTS PREDICTED FROM DISF 

FOR THE 20 AUGUST 1975 CASE 

Key: Model Upwind relocation 

Measured data Downwind relocation 
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FIGURE 24 PREDICTED PATHS OF MAXIMUM CONCENTRATION AT 
600 m ALTITUDE IN THE 10 DECEMBER 1974 CASE 


Shaded area represents region where airborne samples were taken 









MAXIMUM HCI DOSAGE , ppr 



FIGURE 26 COMPARISONS OF GROUND LEVEL HCI MAXIMUM DOSAGE PREDICTIONS 
FOR THE 20 AUGUST 1975 CASE USING OBSERVED CLOUD AS INPUT 
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FIGURE 27 COMPARISONS OF AIRBORNE SAMPLING DATA WITH RESULTS PREDICTED FROM DISF 
USING OBSERVED CLOUD AS INPUT FOR THE 20 AUGUST 1975 CASE 

Key: Model 

Measured data 

Downwind relocation 











KN 0 t=0 


O,# r = o.i , a = 2.00 
A, A r = o.i , a = io.oo 
<£> r = o.oi, a = 2.oo 



Curves are analytical solutions (Eq. (93)). 
a = mass class size; x = time step. 



WITH DIAM< d, ml 


n ii mu 



FIGURE 29 COMPARISON OF COMPUTED AND ANALYTICAL SOLUTIONS 
TO THE COAGULATION EQUATION (EQ. (90)) 


174 


Curves are analytical solutions (Eq. (94)). 
a = mass class size; t = time step. 



WITH DIAM < d, ml 



Curves are drawn through points calculated with the smallest 
time step and mass class sizes at each time, 
a = mass class size; x = time step. 
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FIGURE 31 PARTICLE SIZE AND MASS DISTRIBUTIONS 
FOR THE 20 MAY 1975 CASE 

Dawborn's Alumina Distribution with Debris Having 
1000 Times the Mass of Alumina 








FIGURE 33 PARTICLE SIZE AND MASS DISTRIBUTIONS 
FOR THE 20 MAY 1975 CASE 
Kreautle's Alumina Distribution with Debris 
Having 1000 Times the Mass of Alumina 
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FIGURE 35 MASSES OF ALUMINA AND 
GROUND CLOUD DURING CLOUD 

DEBRIS LOST FROM 
RISE PERIOD 

Mass Ratio of Debris to Alumina = 1000 
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